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Abstract 

The  Generalized  Lambda  Distribution  (GLD)  is  a  four  parameter,  continuous  function 
which  is  flexible  enough  to  approximate  many  different  probability  density  functions  (PDF). 
Like  any  distribution,  the  GLD  can  be  fitted  to  sample  data,  and  can  be  used  to  generate 
similar  samples  for  use  in  computer  simulations  and  models.  The  strengths  of  the  GLD  lie  in 
its  abilities  to  approximate  many  well-known  distributions,  represent  data  when  the  underly¬ 
ing  distribution  is  unknown,  and  simply  and  efficiently  fit  or  generate  random  variates.  The 
method  of  moments  is  presently  the  most  widely-accepted  technique  for  estimating  this  dis¬ 
tribution  from  sample  data.  However,  sample  moments  are  sensitive  to  extreme  observations 
and  are  subject  to  large  sampling  variability  as  the  size  of  the  sample  data  decreases. 

L-moments  are  expectations  of  certain  linear  combinations  of  order  statistics.  They 
can  be  defined  for  any  random  variable  whose  mean  exists.  L-moments  form  the  basis  of  a 
general  theory  that  summarizes  and  describes  probability  distributions  and  observed  data 
samples.  They  can  be  used  to  estimate  parameters  and  quantiles  of  probability  distributions 
as  well.  Their  main  advantage  over  conventional  moments  is  that,  being  linear,  they  suffer 
less  from  the  effects  of  sampling  variability,  and  are  theoretically  more  robust  to  outliers  in 
the  data  and  small  sample  sizes. 

The  process  of  estimating  the  parameters  of  a  probability  distribution  by  matching 
its  L-moments  to  those  corresponding  to  sample  data  is  referred  to  as  the  method  of  L- 
moments.  This  method  has  been  successfully  applied  to  many  well-known  probability  distri¬ 
butions.  In  this  thesis,  the  method  of  L-moments  was  applied  to  estimation  of  the  General¬ 
ized  lambda  distribution.  A  Monte  Carlo  experiment  was  performed  to  compare  the  method 
of  L-moments  to  the  conventional  method  of  moments  and  a  method  which  uses  alternate 
measures  of  symmetry  and  tailweight. 

Experiment  results  showed  that  the  method  of  L-moments  is  better  than  the  con¬ 
ventional  and  alternate  methods  of  moments  for  fitting  distributions  to  sample  data.  The 


xi 


improvement  was  more  noticeable  as  the  skewness  and  kurtosis  of  the  empirical  distribution 
of  sample  data  increased.  There  was  no  change  in  the  advantage  of  L-moments  as  the  sample 
size  changed.  The  experiment  proved  that  L-moments  suffer  less  from  the  effects  of  sampling 
variability,  and  that  they  axe  more  robust  to  outliers  in  the  data,  when  used  to  fit  the  GLD 
to  sample  data. 
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THE  USE  OF  L-MOMENTS  TO  FIT 
THE  GENERALIZED  LAMBDA  DISTRIBUTION 
TO  SAMPLE  DATA 


I.  Introduction 

1.1  Simulation  and  the  Lambda  Distribution 

Industry,  commerce,  government,  and  society  in  general  all  face  problems  which  con¬ 
tinue  to  grow  in  size  and  complexity.  There  is  an  obvious  need  for  timely,  cost-effective  tools 
for  solving  these  problems.  One  such  tool  is  the  computer-aided  simulation  model.  Because 
of  the  rapid  growth  in  the  power  and  abilities  of  small  computers,  computer  simulation  mod¬ 
els  of  dynamic,  complex,  or  unique  situations  are  more  practical  and  easier  to  implement 
than  ever. 

The  problem  with  simulation  is  that  many  applications  require  probabilistic  models  of 
input  variables  based  on  real-world  data.  Data  such  as  the  expected  lifetime  of  a  machine 
part,  the  number  of  troops  needed  to  secure  an  objective,  or  the  radiation  pattern  from  a 
nuclear  explosion  are  needed  by  the  simulation  to  get  results  that  “mean  something.”  How¬ 
ever,  such  information  is  rarely  known  with  certainty.  Instead,  a  planner  may  have  a  bank  of 
accumulated  data  from  previous  tests  and  experiences  which  he  can  use  to  probabilistically 
describe  the  range  of  possible  values  for  each  of  his  variables  of  interest  (Bergevin,  1993:  2). 

The  probability  density  function  (PDF)  summarily  describes  a  continuous  random 
variable.  Many  different  PDF’s,  including  those  of  well-known  distributions  as  the  Normal, 
Exponential,  Gamma,  and  Weibull,  are  available  to  the  modeler.  Each  PDF  has  its  own  dis¬ 
tinct  shape,  or  family  of  shapes,  and  covers  a  specific  range  of  values.  Since  most  commercial 
simulation  packages  provide  a  collection  of  PDF’s,  which  together  cover  the  ga  of  these 
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well-known  random  variables,  the  modeler  must  decide  which  of  the  available  distributions 
best  “fits,”  or  describes,  his  set  of  data. 

An  alternative  to  this  collection  of  specialized  PDF’s  is  a  distribution  with  a  generalized 
density  function  that  is  not  limited  to  a  particular  shape,  but  can  instead  take  on  a 
variety  of  shapes.  Using  such  a  generalized  function  frees  the  modeler  from  choosing  an. 
competing  distributions  and  allows  him  to  possibly  create  an  even  better  fit  to  the  data. 

The  Generalized  Lambda  Distribution  (GLD)  has  such  a  density  function.  The  GLD  is 
a  four  parameter,  continuous  distribution  whose  PDF  is  flexible  enough  to  approximate  many 
different  probability  density  functions.  The  four  parameters,  Aj,  A2,  A3,  and  A4,  together 
define  the  location,  scale,  and  shape  of  the  density  function.  Like  any  distribution,  the 
GLD  can  be  fitted  to  sample  data  and  used  to  generate  similar  samples  for  in  computer 
simulations.  The  strengths  of  the  GLD  are  its  ability  to  approximate  many  well-known 
distributions,  represent  data  when  the  underlying  distribution  is  unknown,  and  simply  and 
efficiently  fit  or  generate  random  variates  (Ramberg  et  ah,  1979:  201). 

The  method  of  moments  is  the  most  widely-accepted  technique  for  estimating  the 
parameters  of  this  distribution  from  sample  data.  Simply  put,  this  method  computes  the  first 
four  sample  moments:  the  mean,  variance,  skewness,  and  kurtosis.  Then,  a  computerized 
search  routine  finds  the  GLD  parameters  that  equate  to  those  moments.  Unfortunately, 
the  higher-order  sample  moments,  the  skewness  and  kurtosis,  are  third-  and  fourth-power 
functions,  respectively,  and  can  vary  greatly  from  sample  to  sample,  especially  when  the 
sample  size  is  small  (Ramberg  et  ah,  1979:  205).  Recent  research  suggests  that  using  Linear 
moments  (L-moments)  in  place  of  conventional  moments  could  eliminate  the  instability  of 
the  GLD  brought  on  by  the  inherent  variability  of  the  higher-order  c<  entional  moments 
(Bergevin,  1993:  4). 

L-moments  are  expectations  of  linear  combinations  of  order  statistics.  They  can  be 
defined  for  any  random  variable  whose  mean  exists.  L-moments  form  the  basis  of  a  theory 
that  summarizes  and  describes  probability  distributions  and  observed  data  samples.  They 
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can  be  used  to  estimate  parameters  and  quantiles  of  probability  distributions  as  well.  Their 
main  advantage  over  conventional  moments  is  that,  being  linear,  they  suffer  less  from  the 
effects  of  sampling  variability,  and  are  theoretically  more  robust  to  outliers  in  the  data 
(Hosking,  1990:  105). 

Estimating  the  parameters  of  a  distribution  by  matching  its  L-moments  to  those  cor¬ 
responding  to  sample  data  is  referred  to  as  the  method  of  L-moments.  This  method  has 
been  successfully  applied  to  many  well-known  probability  distributions,  such  as  the  Normal, 
Gamma,  Uniform,  Exponential,  Laplace,  Logistic,  Bernoulli,  Wakeby  and  Weibull  Distribu¬ 
tions  (Hosking,  1992:  186). 

1.2  Applying  Linear  Moments 

In  theory,  the  method  of  L-moments  could  provide  accurate  estimates  of  the  four 
lambda  parameters  of  otherwise  troublesome  data  samples.  With  accurate  parameter  es¬ 
timates,  the  GLD  should  produce  probability  distribution  functions  for  data  samples  from 
a  wide  variety  of  well-known  distributions.  In  addition,  the  generated  probability  density 
function  (PDF)  should  score  well  on  accepted  goodness-of-fit  tests.  The  end  result  would 
be  a  GLD  that  produces  an  acceptable  fitted  PDF  given  any  data  sample,  including  small 
samples  which  previously  produced  poor  estimates. 

To  date,  researchers  have  proven  the  usefulness  of  the  GLD  in  simulating  other  well- 
known  distributions  (Bergevin,  1993:  3).  The  GLD  has  been  proven  reliable  when  estimating 
distributions  from  large  data  samples  using  the  method  of  moments.  Other  researchers  have 
also  proven  the  usefulness  of  the  method  of  L-moments  in  estimating  the  parameters  of 
well-known  distributions  from  small  data  samples  (Hosking,  1990:  112).  However,  no  one 
has  used  the  method  of  L-moments  to  estimate  the  parameters  of  the  generalized  lambda 
distribution,  particularly  with  small  data  samples,  although  Hosking  (1986)  provides  some 
initial  work  in  this  direction. 


1.3  Overview 


The  purpose  of  thesis  is  to  assess  the  usefulness  of  the  method  of  L-moments  applied 
to  the  GLD,  especially  in  relation  to  other,  more  well-established  methods  for  parameter 
estimation.  One  specific  goal  is  to  determine  whether  any  one  method  is  superior  in  fitting  the 
generalized  lambda  distribution  to  sample  data  when  the  underlying  distribution  is  unknown 
and  the  sample  size  is  relatively  small. 

Chapter  Two  provides  the  background  supporting  this  thesis.  It  presents  the  history 
of  the  GLD,  methods  of  estimating  the  GLD  from  sample  data,  and  an  introduction  to  linear 
moments. 

Chapter  Three  reviews  a  variation  of  the  Lambda  Distribution  proposed  by  Hosking.  It 
examines  and  reconciles  the  differences  between  Hosking’s  five-parameter  Lambda  Distribu¬ 
tion  and  the  four-parameter  Generalized  Lambda  Distribution  of  Ramberg,  Schmeiser,  and 
Tukey.  A  resolution  of  the  differences  in  the  L-moment  equations  used  for  the  two  different 
distributions  is  presented.  Finally,  it  compares  the  estimation  techniques  of  linear  moments 
developed  by  Hosking  to  two  suggested  by  Mykytka  and  Bergevin. 

Chapter  Four  compares  the  search  routines  used  to  solve  the  simultaneous  equations 
which  determine  the  lambda  parameters  from  the  sample  L-moments.  Hosking  uses  a 
Newton-Raphson  iteration  to  find  his  five  lambda  parameters,  while  Mykytka  and  Ram¬ 
berg  use  an  adaptation  of  Powell’s  algorithm  to  find  the  four  parameters  of  the  GLD.  The 
comparison  examines  it  the  accuracy  of  the  algorithms,  the  computational  speed  of  the 
routines,  and  the  robustness  of  the  functions,  for  each  of  the  two  methods. 

Chapter  Five  discusses  the  design  of  a  Monte  Carlo  experiment  used  to  compare  the 
effects  of  sampling  variability  on  the  method  of  linear  moments,  the  method  of  conventional 
moments,  and  the  method  of  alternate  moments.  The  effectiveness  of  the  three  methods  was 
evaluated  based  on  visual  comparisons  of  PDF’s  and  Kolmogorov-Smirnov  goodness-of-fit 
statistics.  The  experiment  compared  the  three  methods  over  three  different  distribution 
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shapes:  Normal,  Gamma,  and  Exponential;  and  three  different  sample  sizes:  n  =  25,  n  =  50, 
and  n  =  100. 

Chapter  Six  compares  test  results  of  the  method  of  L-moments  with  those  of  the  con¬ 
ventional  method  of  moments  and  an  alternate  method  proposed  by  Mykytka.  It  compares 
the  goodness-of-fit  using  both  Kolmogorov- Smirnov  test  statistics  and  visual  plots.  It  also 
examines  the  method’s  stability  with  respect  to  small  data  samples  by  experimenting  with 
sample  sizes  of  100,  50,  and  25  elements.  Finally,  it  examines  the  method’s  flexibility  with 
respect  to  different  distribution  shapes  by  testing  samples  from  the  GLD’s  approximations 
to  the  Normal,  Gamma,  and  Exponential  distributions. 

Chapter  Seven  discusses  the  conclusions  drawn  from  the  thesis,  and  makes  recommen¬ 
dations  for  further  research  studies. 
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II.  Discussion  of  the  Literature 

2.1  Chronological  Organization 

This  chapter  chronologically  reviews  major  research  milestones  for  both  the  Gener¬ 
alized  Lambda  Distribution  and  Linear  Moments.  It  begins  with  the  history  and  recent 
development  of  the  lambda  distribution.  Although  the  L-moments  were  developed  at  the 
same  time  as  the  GLD,  they  were  only  recently  applied  to  the  generalized  lambda  distribu¬ 
tion.  It  is  from  this  point  of  discovery  that  the  discussion  of  L-moments  begins. 

2.2  History  of  the  Lambda  Distribution 

Although  most  continuous  probability  distributions  are  defined  in  terms  of  their  density 
functions,  f(x),  or  cumulative  distribution  functions,  F(x),  it  is  valid  to  define  a  distribution 
by  its  percentile  function,  if  that  percentile  function  exists.  The  percentile  function  is  simply 
the  inverse  of  the  distribution  function,  i.e.,  R(p)  —  or  equivalently,  p  =  F(x).  The 

percentile  function  determines  the  value,  x  =  M(p),  such  that  the  probability  that  a  random 
variable  having  this  distribution  takes  on  a  value  less  than  x  is  p  (Mykytka,  1979:  362).  This 
relationship  is  depicted  in  Figure  2.1. 

The  ability  to  express  a  random  variable  in  terms  of  its  percentile  function  is  quite 
useful  in  Monte  Carlo  simulation  studies.  In  particular,  it  is  well  known  that  if  R  is  the 
percentile  function  of  a  continuous  probability  distribution  and  if  the  random  variable  U 
is  uniformly  distributed  on  the  (0,1)  interval,  then  the  transformation  X  =  R{U)  yields 
a  continuous  random  variable  with  percentile  function  R.  Thus,  since  sources  of  uniform 
(0,1)  pseudo-random  variates  are  commonly  available,  this  transformation  yields  a  simple 
method  for  generating  pseudo-random  variates  from  distributions  whose  percentile  functions 
are  known  and  are  computationally  tractable  (Ramberg  et  al.,  1979:  202). 
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Figure  2.1  The  relationship  between  the  density  function  and  the  percentile  function 

Tukey  (1960)  created  a  function,  which  he  called  the  lambda  function,  with  this  form. 
Tukey’s  function,  which  is  valid  for  all  non- zero  A,  can  be  written  as 

—  ( 1  —  p )* 

R(p)  =  - — (2-1) 

In  1976,  Ramberg  et  al.  (1979),  generalized  Tukey’s  single- parameter  Lambda  distri¬ 
bution  and  named  it  the  Generalized  Lambda  Distribution.  They  recognized  its  strengths, 
by  stating  that  “important  applications  of  the  distribution  include  the  modeling  and  sub¬ 
sequent  generation  of  random  variates  for  simulation  studies  and  Monte  Carlo  sampling 
studies  of  the  robustness  of  statistical  procedures”  (Ramberg  et  al.,  1979:  201).  Mykytka 
(1978)  further  noted  that  “defined  by  its  percentile  function,  the  GLD  provides  a  simple  and 
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computationally  efficient  method  for  generating  random  variates  having  a  wide  variety  of 
distributions.” 

Ramberg  et  al.  generalized  Tukey’s  function  given  by  Equation  2.1  to  a  four-parameter 
form  that  could  be  used  to  approximate  a  number  of  well-known  symmetric  and  asymmetric 
distributions.  Their  distribution  is  defined  by  the  percentile  function 


R(p)  —  Ai  + 


(i  -  pY 


-,  0  <  p  <  1. 


Ramberg  et  al.  note  that  a  close  approximation  to  the  standard  normal  distribution 
results  when  Ai  =  0,  A2  =  0.1975,  and  A3  =  A4  =  0.1350  (Ramberg  et  al.,  1979:  203).  The 
distribution  defined  by  Equation  2.2  is  referred  to  as  the  Generalized  Lambda  Distribution. 
The  GLD  has  also  been  referred  to  as  the  Ramberg-Schmeiser-Tukey  (RST)  distribution. 
The  parameters  At  and  A2  are  location  and  scale  parameters,  respectively,  while  A3  and  A4 
are  shape  parameters  that  jointly  determine  the  skewness  and  kurtosis  of  the  GLD.  When 
A3  =  A4,  the  resulting  density  is  symmetric. 

Using  the  fact  that  F(x )  =  p  and  x  =  R(p),  we  can  find  the  density  function  corre¬ 
sponding  to  Equation  2.2  by  noting 


f(x)  = 


dF(x ) 


dp 

dm 


dR{p) 

dp 


which  yields 


d_m 

dp 


A 3pX3_1  +  A4(l  -  p)x*~ 1 


0  <  p  <  1. 


The  cumulative  distribution  function  of  the  GLD  does  not,  in  general,  exist  in  a  simple 
closed  form,  but  this  is  not  a  cause  for  concern  since  it  is  also  true  of  the  Normal  distribution, 
whose  percentiles  are  more  difficult  to  compute.  For  the  GLD,  it  is  simple  to  obtain  plots  of 
the  distribution  function  bv  plotting  p  on  the  y-axis  versus  R(p)  on  the  x-axis.  Similarly,  a 


plot  of  the  density  function  is  obtained  by  plotting  f[R(p)]  on  the  y-axis  against  R(p)  on  the 
x-axis,  for  p  ranging  from  zero  to  one.  FORTRAN  programs  that  compute  R(p)  and  f[R(p)] 
for  specified  lambda  values  are  given  in  Mykytka  (1978:  82-84)  and  in  the  appendices  of  this 
study. 

Ramberg  and  Schmeiser  developed  the  following  expressions  for  the  mean,  variance, 
skewness,  and  kurtosis  of  the  GLD: 


M(Al  >  A2,  A3,  A4)  —  +  — , 

A2 


<7 


2 


<72(A2,  A3,  A4)  — 


B  —  A2 

A22 


«3 


«3(A3,  A4)  =  sign( A2)  • 


C  -  3 AB  +  2A3 
{B-A2)i  ’ 


a4 


o4(A3,  A4) 


D  -  AAC  +  6 A2B  -  3 A4 
(B  —  A2)2 


where 


1  _ 1 

1  +  A3  1  +  A4 


(2.4) 


(2.5) 

(2.6) 

(2.7) 


(2.8) 


B  —  - - — - 2/7(1  A3, 1  +  A4)  -f 


1  -f  2A3 


1  +  2A4  ’ 


C  —  - — - rr  3/9(1  +  2A3, 1  +  A4)  +  3/9(1  +  A3, 1  +  2A4) 


1  +  3A3 


1  -f  3A4 


(2.9) 


(2.10) 


D  = 


- — ■ — — - 4,9(1  +  3A3, 1  +  A4)  +  6/9(1  2A3, 1  +  2A4) 

1  +  4A3 

—4/7(1  +  A3, 1  +  3A4)  +  .  , 

1  +  4  A., 


(2.11) 
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Figure  2.2  Characterization  of  Various  Distributions  by  Skewness  and  Kurtosis 

and  sign( A2)  will  be  either  1  or  —1,  depending  on  whether  the  A2  parameter  is  positive  or 
negative. 

As  the  notation  indicates,  the  skewness  and  kurtosis  are  functions  of  A3  and  A4  alone. 
The  skewness,  03,  is  a  function  of  only  A3  and  A4  since,  as  long  as  A3,  A4  >  0,  the  sign  of  A2  is 
always  the  same  as  the  sign  of  both  A3  and  A4  (Mykytka,  1978:  20).  The  variance,  however, 
also  depends  on  the  shape  parameter  A2  and  the  mean  depends  on  all  four  parameters.  The 
shaded  region  of  Figure  2.2  shows  the  different  combinations  of  skewness  and  kurtosis  that 
can  be  obtained  from  the  GLD. 

The  technique  known  as  the  method  of  moments  is  the  usual  means  of  selecting  the 
values  of  the  GLD  parameters.  The  method  of  moments  begins  by  estimating  the  mean, 
variance,  skewness,  and  kurtosis  of  a  sample.  Skewness  and  kurtosis  are  the  standardized 
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third  and  fourth  moments  as  defined  by 


<*3 


4 


—  E(x~»f 

<73  ’ 

_  E(x-n)4 
a4 


(2.12) 

(2.13) 


Since  the  values  for  fi,  a2,  a3,  and  a4  are  unknown,  they  must  be  estimated  by  the 
sample  statistics  X,  S'2,  a3,  and  d4,  respectively.  These  sample  moments  are  defined  by 


X 

=  iai^ 

(2.14) 

s' 2  = 

=  ^=l(Xi-X)\ 

(2.15) 

-  £E"q(S,-X)3 

S' 3  ’ 

(2.16) 

d?4 

_  iET-i  (*-*)« 
s' 4 

(2.17) 

Once  the  first  four  sample  moments  are  computed,  they  arc  applied  to  Equations  2.4 
through  2.11  to  solve  for  the  lambda  parameters.  Further,  by  choosing  the  four  parameter 
values  appropriately,  a  wide  range  of  distributions  can  be  duplicated,  as  indicated  by  the 
shaded  portion  of  Figure  2.2 

The  biggest  problem  for  the  OLD,  though,  is  the  presence  of  variables  raised  to  the  third 
and  fourth  power,  respectively,  in  the  Equations  for  the  third  and  fourth  sample  moments. 
Ramberg  et  al.  said  that  “sample  moments  are  sensitive  to  extreme  observations,  such  as 
values  more  than  two  standard  deviations  from  the  mean,  and  sampling  variability  of  the 
third  and  fourth  moments  can  be  large.”  (Ramberg  et  al.,  1979:  205).  This  can,  in  turn, 
result  in  a  poor,  or  even  unacceptable,  GLD  fit.  Mykytka  recommends  that  the  method  of 
moments  be  used  to  fit  data  only  when  the  sample  size  exceeds  150.  (Mykytka.  1978:  00). 


He  also  suggests  that  alternatives  to  the  method  of  moments,  like  non-linear  least  squares 
or  using  four  specified  percentiles  to  solve  for  the  parameters  simultaneously,  be  explored 
(Mykytka,  1978:  80). 


2.3  Estimating  Distributions  from  Data 

Mykytka  and  Ramberg  acted  upon  their  recommendations  and  experimented  with  the 
GLD  using  an  alternative  to  the  method  of  moments.  They  developed  a  new  procedure  using 
two  statistics  that  are  functions  of  order  statistics  and  compared  it  to  the  method  of  moments 
(Mykytka,  1979:  361).  The  two  statistics  employed  in  the  alternative  procedure  provide 
measures  of  the  symmetry  and  tail  weight  of  a  distribution.  Hence,  they  are  analogous 
to  the  third  and  fourth  moments,  i.e.,  skewness  and  knrtosis.  The  two  procedures  were 
compared  by  a  Monte  Carlo  study  which  examined  the  effects  of  sampling  variability  on 
them  (Mykytka,  1979:  362). 

For  this  effort,  Mykytka  and  Ramberg  chose  two  statistics  developed  by  Hogg,  Fisher, 
and  Randles  (1975)  that  are  functions  of  linear  combinations  of  the  order  statistics.  To 
define  these  statistics,  let  Y\  <  V-2  <  ...  <  V„  be  the  order  statistics  of  a  random  sample  of 
size  n  from  a  continuous  distribution.  The  skewness  indicator,  Q3 ,  is  then  defined  by 


0  (0.05)— M(0.50) 
“  M (0.50)  — L(0.05) : 


(2.18) 


where  C/ (0.05),  MfO.50),  and  Z(0.05)  are,  respectively,  the  averages  of  the  largest  (upper¬ 
most)  five  percent  of  the  Y's,  the  middle  fifty  percent  of  the  Y's,  and  the  smallest  (lowest) 
five  percent  of  the  Y's.  Similarly,  Q4,  the  taihveight  descriptor,  is  defined  as 


n  _  0 (0.05)— L(0.05) 
—  t/(0.50)-L(0.50)* 


(2.19) 
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These,  combined  the  first  two  conventional  moments,  the  sample  mean  and  variance,  com¬ 
prised  the  alternative  method  of  estimation  (Mykytka  and  Ramberg,  1979:  364).  Mykytka’s 
alternative  method  proved  more  accurate  at  estimating  the  GLD  parameters  of  sample  data 
than  the  conventional  method  of  moments.  From  their  Monte  Carlo  study  they  concluded 
that  “the  curve  shapes  of  the  densities  estimated  by  the  method  of  moments  are  more  varied 
than  the  shapes  of  the  densities  estimated  by  the  [alternate]  method.”  This  reduced  vari¬ 
ation  in  estimated  densities  suggests  that  sampling  variability  affects  the  alternate  method 
less  than  it  affects  the  method  of  moments.  This  alternative  effectively  expanded  the  range 
of  allowable  sample  sizes  because  the  alternative  method  more  consistently  estimated  pa¬ 
rameters  for  samples  as  small  as  fifty  elements,  as  opposed  to  the  method  of  moments,  which 
was  inconsistent  even  for  samples  of  100  elements  (Mykytka,  1979:  366). 

Ozturk  and  Dale  (1985)  developed  and  investigated  the  method  of  nonlinear  least 
squares  estimation,  based  on  the  suggestion  of  Ramberg  et  al.  (1979).  They  used  least 
squares  estimation  procedures  to  estimate  the  parameters  of  the  GLD  by  minimizing  the 
squared  distance  between  the  theoretical  and  empirical  percentile  functions.  Although  the 
least  squares  procedure  offered  no  major  computational  advantage  over  the  method  of  mo¬ 
ments,  it  was  easier  to  compute  the  objective  function  and  the  first  and  second  lambda 
parameters  (Ozturk  and  Dale,  1985:  81-82).  Ozturk  and  Dale  said  that,  “to  compare  the 
accuracy  and  precision  of  the  fits  obtained  by  the  different  methods,  we  performed  a  Monte 
Carlo  simulation  study  similar  to  that  made  by  Mykytka  and  Ramberg  (1979).  The  results 
showed  no  appreciable  differences,  as  long  as  the  parameters  could  be  estimated.”  (Ozturk 
and  Dale,  1985:  83).  The  method  of  moments  occasionally  produces  sample  values  of  skew¬ 
ness  and  kurtosis  for  which  the  corresponding  A3  and  A.t  parameters  do  not  exist.  Ozturk 
and  Dale  (1985)  concluded  that  the  least  squares  estimation  procedure  does  not  suffer  from 
this  problem. 

Cheng  (1985)  developed  an  alternate  procedure  using  the  sample  median  and  interquar¬ 
tile  range  to  select  the  parameter  values.  She  compared  it  to  the  method  of  moments  by 
means  of  a  Monte  Carlo  simulation  study.  Cheng  concluded  that  the  alternate  procedure 
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appeared  to  be  no  better,  and  in  some  cases,  slightly  worse,  than  the  method  of  moments. 
She  also  extended  the  concept  to  a  procedure  to  select  lambda  parameter  values  which  match 
specified  percentiles  of  the  distribution,  but  concluded  that  it  was  no  more  promising  than 
her  first  procedure.  Cheng  used  the  chi-square  goodness-of-fit  test  to  compare  the  different 
methods,  but  had  difficulty  determining  the  actual  degrees  of  freedom  for  the  tests.  She 
hypothesized  that  the  problem  was  caused  by  a  “large”  number  of  parameters  in  the  GLD. 
In  this  case,  the  test  statistic  may  not  possess  an  approximate  chi-square  distribution,  or 
the  method  of  estimating  the  number  of  degrees  of  freedom  is  inappropriate  for  this  distri¬ 
bution.  The  fact  that  neither  the  method  of  moments  nor  the  alternate  method  produce 
maximum-likelihood  estimators  also  caused  difficulty  with  the  chi-square  test  (Cheng,  1985: 
52-54). 

Chou  (1988)  developed  the  method  of  maximum  likelihood  for  determining  the  four 
lambda  parameters  of  the  GLD.  She  compared  her  method  to  the  method  of  moments 
to  determine  which  method  produces  better  fits  and  better  parameter  estimates.  Chou 
conducted  a  Monte  Carlo  experiment  to  perform  the  comparison.  Her  experimental  results 
showed  that  the  method  of  maximum  likelihood  appears  to  produce  slightly  better  fits  and 
better  parameter  estimates  than  the  method  of  moments.  Unfortunately,  an  in-depth  study 
was  restricted  by  the  computational  demands  of  the  MLE  method.  Chou  reports  that  the 
computer  CPU  time  used  by  the  MLE  method  is  about  one  hundred  times  more  than  that 
used  by  the  method  of  moments,  and  the  CPU  time  needed  grows  exponentially  with  the 
number  of  observations  in  the  sample.  She,  therefore,  recommended  the  continued  use  of 
the  method  of  moments  (Chou,  1988:  70-71). 

2.4  Computer  Software  for  the  Lambda  Distribution 

By  1991,  the  computer  revolution  had  created  an  environment  where  work  requiring 
a  mainframe  computer  in  1979  could  be  done  on  a  home  computer.  The  state  of  the  art 
in  computer  languages  had  progressed  from  FORTRAN  77  to  object-oriented  languages  like 

C+  +  . 
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Chung-Lung  Hsu  wrote  a  computer  program  in  C++  in  1991,  as  a  Master’s  thesis  effort, 
which  implemented  methods  for  finding  values  for  the  lambda  parameters.  The  methods 
computed  parameters  either  from  sample  data  or  from  various  statistics,  and  used  the  lambda 
parameters  to  plot  the  curve  of  the  probability  distribution  (Hsu,  1991:  2).  When  given 
sample  data  or  statistics  of  the  data,  the  program  used  these  techniques  to  compute  the 
lambda  parameters:  the  method  of  moments;  the  method  of  moments  with  boundary  points 
known;  the  method  of  percentiles;  the  replacement  of  mean  and  variance  with  the  median 
and  interquartile  range;  the  replacement  of  skewness  and  kurtosis  with  alternate  measures  of 
symmetry  and  tailweight  (Qz,Qa)‘,  and  the  method  of  least  squares  estimation  (Hsu,  1991: 
13-34). 

All  of  these  methods  were  suggested  by  Mykytka  (Mykytka,  1978:  80-81)  or  Ramberg 
(Ramberg  et  al.,  1979:  205).  In  fact,  many  of  Hsu’s  C++  algorithms  were  adaptations  of 
Mykytka’s  FORTRAN  algorithms  written  in  1979.  Hsu  states  that  “all  known  methods  of 
determining  the  lambda  parameters  have  been  introduced.”  (Hsu,  1991:  29).  Despite  the 
number  of  alternatives  attempted,  Hsu  concluded  that  none  was  a  significant  improvement 
over  the  original  method  of  moments  first  proposed  by  Ramberg  (Hsu,  1991:  29-33). 

2.5  Analysis  of  the  Generalized  Lambda  Distribution 

Bergevin  (1993)  analyzed  the  generalized  lambda  distribution  intent  on  expanding  the 
range  of  permissible  skewness  and  kurtosis  values  used  to  compute  lambda  parameters.  An 
increase  in  the  parameter  range  increases  the  number  and  types  of  probability  distributions 
that  the  GLD  can  approximate  (Bergevin,  1993:  16).  While  he  was  not  able  to  expand 
the  coverage  range  of  the  GLD  over  all  possible  values  of  skewness  and  kurtosis,  Bergevin 
successfully  reduced  the  previously-uncovered  range  of  values  by  approximately  fifty  percent 
(Bergevin,  1993:  53). 

Because  of  its  ability  to  “match  moments,”  the  GLD  can  mimic  the  behavior  of  most 
common  PDF’s  by  setting  the  four  parameters  to  appropriate  values.  Figure  2.2  shows  the 
ranges  of  these  PDF's  in  terms  of  their  measures  of  skewness  and  kurtosis.  Some  distri- 
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butions,  such  as  the  Uniform,  Normal  and  Exponential,  are  represented  by  single  points: 
others,  such  as  the  Student’s  t,  Log-Normal  and  Gamma,  are  represented  by  curves;  the 
three  types  of  beta  distributions  are  represented  by  regions  of  values.  The  top  right-hand 
region,  denoted  the  “Impossible  Area,”  contains  skewness-kurtosis  combinations  that  will 
never  be  exhibited  by  any  PDF.  The  GLD,  at  present,  can  “mimic”  those  P£)F’s  found 
within  the  shaded  regions  of  Figure  2.2.  Obviously,  all  but  some  distributions  with  moments 
similar  to  those  of  the  U-shaped  beta  distributions  can  currently  be  modeled  using  the  GLD 
(Bergevin,  1993:  53). 

The  portion  of  his  research  that  is  critical  to  this  study  is  the  theory  of  L-moments. 
During  an  on-going  literature  search  in  the  latter  stages  of  his  thesis  effort,  Bergevin  un¬ 
covered  publications  of  Hosking’s  current  research  of  a  new  class  of  linear  statistics,  which 
Hosking  named  L-moments.  While  there  was  insufficient  time  to  thoroughly  develop  an  ap¬ 
plication  of  L-moments  to  the  GLD,  Bergevin  conducted  a  promising  initial  analysis  of  the 
idea  and  derived  the  equations  that  express  the  theoretical  L-moments  as  functions  of  the 
four  lambda  parameters  (Bergevin,  1993:  59-60). 

2.6  Derivation  of  the  GLD  from  L-Moments 

Hosking  defines  the  first  four  L-moments  via 


A, 

=  /  R(p)dp, 

Jo 

(2.20) 

A  2 

=  f  R(p )  •  (2p  —  l)dp, 

Jo 

(2.21) 

A3 

=  [  R{P)  ■  (6P2  -  6p  +  1  )dp, 

Jo 

(2.22) 

a4 

=  C  R{p )  ■  (20p3  -  30p2  +  12p  +  1  )dp, 

Jo 

(2.23) 

where  R(p)  is  simply  the  percentile  function. 

As  with  the  commonly-used  measures  of  skewness  and  kurtosis,  Hosking  (1986:  T) 
chose  to  define  the  two  higher  order  moments  as  dimensionless  ratios  relative  to  the  second 
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order  moment 


_  A a 
A2  * 


U 


This  research  conforms  to  that  convention. 

The  four  L-moments  have  roles  similar  to  those  of  conventional  moments.  Since  the 
first  L-moment  is  simply  the  expected  value  of  the  distribution,  it  is  identical  to  the  mean 
of  the  distribution.  Also,  a  symmetric  distribution  will  have  r3  =  0,  just  as  a3  =  0  for 
symmetric  distributions  using  conventional  moments.  According  to  Hosking,  however,  the 
estimators  of  r3  and  r4  are  more  stable  measures  than  c?3  and  a4,  and  are  better  estimates 
of  symmetry  and  tailweight  of  the  underlying  theoretical  PDF  of  a  distribution. 

For  the  GLD,  the  L-moment  Equations  become 


Ai  = 


^1  + 


_ _ Aj-Ag _ . 

A*(A3+1)(A4+1)» 


(2.24) 


^3(^4+ l)(A4-f-2)+ A4  (A3-MHA3-+  2) 
•^2(^3  +  l)(^3  +  2)(A4*fl)(A4+2)  * 


(2.25) 


(A%— A3)(A4+l)(A4-f  2)(A4+3)  — (A^  —  A4)(A3  +  l)(A3+2)(A3+3) 
(A3+3)(A4+3)[A3(A4+l)(A4+2)+A4(A3+l)(A3+2)] 


(2.26) 


_  (Ag-3A3+2A3)(A4+l)(A4+2)(A4+3)(A4+4)+(A;-3Aj+2A4)(A3+l)(A3+2)(A3+3)(A3-H)  (e)  97>. 

^  (A3+3)(A3+4)(A4+3)(A4+4)[A3(A4  +  l)(A4+2)+A4(A3  +  l)(A3+2)]  v  •“  / 

Bergevin’s  thesis  contains  a  complete  derivation  of  these  equations  (1993:  appendix 
A).  Although  Equations  2.24  through  2.27  look  complicated,  they  are  all  composed  of  poly¬ 
nomials.  Unlike  the  GLD  Equations  (Equations  2.4  through  2.11)  based  on  conventional 
moments,  these  have  no  beta  functions.  Since  computerized  solution  algorithms  for  systems 
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of  polynomial  equations  are  fairly  common  and  easily  adaptable,  it  is  easier  to  determine 
the  values  of  the  GLD  parameters  using  L-moments  instead  of  conventional  moments. 

2.7  (The  Problem  of)  Estimating  L-Moments 

Bergevin’s  review  of  the  available  literature  did  not  reveal  how  to  estimate  the  L- 
moments  from  sample  data.  Thus,  he  hypothesized  two  ways  of  doing  so. 

First,  he  noted  that  Hosking  provides  an  alternate  definition  of  the  L-moments  based 
on  the  order  statistics: 


Ax  = 

E(X), 

(2.28) 

A2  = 

\E{X2:2  ~ 

(2.29) 

A3  = 

\E(X 3*  -  2X2.3  +  *1:3), 

(2.30) 

a4  = 

\E{X 4,4  —  3  A  3,4  +  3X2:4  —  A1:4), 

(2.31) 

where  Xa:b  denotes  the  ath  order  statistic  in  a  sample  of  size  b.  Then  he  suggested  that 
one  could  define  estimators  of  these  based  on  the  average  of  the  results  of  every  possible 
subset  of  the  proper  size  within  the  data  set.  For  example,  in  calculating  the  value  A2  for 
a  sample  of  size  n  =  10,  one  could  estimate  the  value  of  A2  by  computing  one-half  of  the 
difference  between  the  largest  and  smallest  observations  within  all  ^2^  possible  subsets  of 
size  2  selected  from  the  10  observations  and  then  averaging  the  results. 

His  second  suggestion  was  to  derive  the  empirical  CDF  of  the  data  set,  and  create 
finite  summations  to  approximate  Equations  2.20  through  2.23.  This  seemed  to  be  the  more 
reasonable  approach,  since  the  work  required  is  less  than  that  involved  in  calculating  each 
possible  subset  of  sizes  2,  3,  and  4. 

If  an  effective  means  of  calculating  sample  L-moments  can  be  found,  then  Equations 
2.24  through  2.27  can  replace  the  more  complicated  GLD  moment  equations  to  estimate 
appropriate  values  for  the  A,.  Of  course,  the  quality  of  these  two  possible  estimation  proce- 
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dures  is  yet  unproven:  no  one  knows  whether  either  method  produces  unbiased  estimators 
or  whether  they  would  provide  a  good  fit  to  the  data  sample  or  population  (Bergevin,  1993: 
60-61).  The  primary  objective  of  this  thesis  is  to  assess  the  latter. 

2.8  Another  Estimation  Method 

In  his  1986  report  which  introduced  the  notion  of  L-moments,  Hosking  proves  mathe¬ 
matically  and  by  example  that  sample  L-moments,  being  linear  functions  of  the  data,  suffer 
less  from  the  effects  of  sampling  variability  than  do  conventional  sample  moments.  They  are 
more  robust  to  outliers  in  the  data,  and  they  enable  more  secure  inferences  to  be  made  from 
small  samples  about  an  underlying  probability  distribution.  He  claims  that  L-moments  fre¬ 
quently  yield  better  parameter  estimates  than  conventional  moments  (Hosking,  1986:  1).  In 
some  cases,  L-moments  can  specify  a  distribution  even  if  some  of  its  conventional  moments 
do  not  exist.  In  fact,  all  linear  moments  exist  if  and  only  if  the  mean  exists  (Hosking,  1986: 
7).  Hosking  concluded  “that  the  main  features  of  a  probability  distribution  should  be  well 
summarized  by  the  following  four  measures:  the  mean  or  L- location,  the  L-scale,  A2;  the 
L-skewness,  73;  and  the  L-kurtosis,  t4  (Hosking,  1986:  9). 

Hosking  further  proposes  a  method  of  estimating  L-moments  from  a  sample,  and  of 
calculating  the  parameters  of  the  GLD  from  those  sample  L-moments.  These  aspects  of  the 
method  of  L-moments  are  discussed  at  length  in  Chapter  Three. 
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III.  Hosking’s  Lambda  Distribution  and  L-Moments 

3.1  Comparing  Ramberg’s  and  Hosking’s  Lambda  Functions 

Hosking  developed  the  L-moment  equations  for  sixteen  specific  distributions,  including 
a  form  of  the  Generalized  Lambda  Distribution.  The  five-parameter  form  of  the  GLD  he 
chose  is  given  by  (Hosking,  1986:  83)  as 

x(F)  =  e  +  o^-7(l  -F)6,  0  <  F  <  1.  (3.1) 


Comparing  Hosking’s  notation  with  that  more  commonly  used  by  Ramberg  et  al.  for  the 
GLD,  and  assuming  that  a  =  7  as  Hosking  does,  reveals  that  F  =  p,  x(F)  —  R(p),  £  =  A1} 
a  =  7  =  yd  =  A3,  8  =  A4.  Substituting  Ramberg’s  notation  in  Equation  3.1  results  in 


R(p)  =  Ai  + 


PX3  ~  (1  -  p)Xt 
A2 


0  <  p  <  I, 


(3.2) 


the  usual  percentile  function  of  the  four- parameter  GLD  (Bergevin,  1993:  7).  Hosking  later 
altered  his  quantile  function  to  ensure  its  continuity  as  /?  and  8  approach  zero.  The  altered 
function,  shown  in  Equation  3.3,  and  its  resulting  L-moment  equations  are  suitable  for  a 
Newton-Raphson  search  algorithm  (Hosking,  1993:  1), 


x(F)  =  i  + 


F0~1  7  '  (1  —  (1  —  F)6) 

0  <5 


0  <  F  <  1. 


(3.3) 


This  form  of  the  Lambda  Distribution  is  implemented  in  Hosking’s  Lambda  Parameter  Es¬ 
timation  FORTRAN  subroutine  found  in  Appendix  N. 


3.2  Comparing  Bergevin’s  and  Hosking’s  L-Moment  Equations 

Using  the  definitions  of  the  first  four  L-moments  given  by  Equations  2.20  through 
2.23,  Bergevin  derived  the  equations  listed  below  for  the  first  four  L-moments  of  the  GLD 
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(Bergevin,  1993:  59): 


.  _  x  , _ A4  ~  A3 

^2(-^3  +  1)(A4  +  1) 

-^3(^4  +  1)(A4  +  2)  +  A4(A3  +  1)(A3  +  2) 
A2(  A3  +  1 )( A3  +  2)(A4  -f  1)(A4  +  2) 


(3.4) 

(3.5) 


A  (A^-A3)(A«-H)(A4+2)(A4+3)-(A;-A4)(A3  +  l)(A3+2)(A3+3) 

3  A2(A3+1)(A3+2)(A3+3)(A4  +  1)(A4+2)(A,  .3) 


(3.6) 


(A3-3A3+2A3)(A4+l)(A4+2)(A4+3)(A4+4)+(A?-3A?+2A4)(A3+l)(A3+2)(A3+3)(A3+4) 
^2(A3+l)(A3+2)(A3+3)(A3+4)(A4  +  l)(A4+2)(A4+3)(A4+4) 


(3.7) 


The  equations  for  the  first  four  L-moments  which  Hosking  derived  from  his  five-parameter 
form  of  the  GLD  are  (Hosking,  1986:  83) 


1 

V— » 

+  -! 

(3.8) 

\  -  al 3  1 

7  6 

(3.9) 

2  (i+m+P) 

(1  +  £)(2  +  sy 

aP(fl  -  1) 

78(8  —  1) 

(3.10) 

(l  +  j8)(2  +  /?)(3  +  0) 

(l  +  $)(2  +  tf)(3  +  «)’ 

0/3(0- l)(0-2)  ■yS(8-l)(S-2) 

(l+0)(2+0)(3+0)(4+0)  f  (n-5)(2+5)(3+6)(4+5)- 


(3.11) 


There  are  apparent  similarities  between  the  L-moment  equations  for  these  forms  of  the  GLD. 
Appendix  A  contains  detailed  proofs  of  the  equivalence  of  each  of  the  four  pairs  of  equations. 


3.3  Comparing  Sample  L-Moments 

Nowhere  in  his  various  papers  and  reports  on  L-moments  does  Hosking  thoroughly 
discuss  estimating  L-moments  from  sample  data.  In  one  of  his  most  specific  statements, 
he  notes  that  the  sample  L-moments  can  be  defined  as  U-statistics  [a  concept  originally 


introduced  by  Hoeffding  (1948)]  as  follows: 


Li  =  n-'^xi,  (3.12) 

i 

L2  =  1(2)  (3.1.3) 

*  v  '  »  3 

U  ~  0(3)  YH2^2(X>  ~  ^X3  +  X*)>  (3-14) 

6  V  '  i>  j>  k 

=  I  (4)  5ZSSE(zi“3xi  +  3xik-*/).  (3.15) 

4  W  i>  j>  k>  l 


The  equation  for  the  first  sample  L-moment  is  simply  the  equation  for  the  sample  mean.  The 
second  equation  implicitly  computes  the  average  of  the  differences  between  the  largest  and 
smallest  observations  within  all  subsets  of  size  2  from  the  sample.  The  third  and  fourth  equa¬ 
tions  similarly  turn  compute  the  average  values  of  certain  functions  of  the  order  statistics 
within  all  subsets  of  sizes  3  and  4,  respectively.  These  four  equations  turn  out  to  imple¬ 
ment  Bergevin’s  suggestion  for  estimating  the  L-moments  based  on  such  order  statistics,  as 
discussed  in  Chapter  Two. 

It  is  not  clear,  however,  that  this  is  the  procedure  actually  implemented  by  Hosking 
since  his  FORTRAN  L-moment  subroutines  (obtained  from  the  Internet  Statistics  Library 
at  Carnegie-Mellon  University  and  listed  in  Appendix  M)  appear  to  have  a  very  different 
structure.  In  fact,  the  first  step  in  Hosking’s  algorithm  computes  what  he  calls  a  “plotting 
position”  which  is  intended  to  compensate  for  the  fact  that  the  finite  data  sample  comes 
from  an  infinite  population  and  allows  for  distributions  with  infinite  tails  (Hosking,  1986: 
16).  This  is  followed  by  a  number  of  manipulations  that  appear  to  be  performed  for  the 
sake  of  computational  efficiency.  Unfortunately,  these  obscure  the  exact  procedure  being 
implemented.  It  does  appear,  however,  that  this  approach  is  tantamount  to  the  computation 
of  the  L-moments  of  the  empirical  CDF,  as  also  suggested  by  Bergevin.  Isolated  comments 
made  by  Hosking  suggest  that  these  two  approaches  are  equivalent. 
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IV.  Lambda  Parameter  Search  Routines 

4-1  Hosking’s  Lambda  Quantile  Function 

Recall  that,  as  discussed  in  Chapter  Three,  Hosking  revised  the  quantile  function  of 
the  five-parameter  version  of  the  GLD  to  have  the  form  shown  in  Equation  3.3.  Prior  to 
this  revision,  the  denominators  of  the  third  and  fourth  L-moments,  r3  and  r4,  were  zero 
when  either  0  ov  8  was  zero,  creating  a  discontinuity.  In  the  revised  form,  since  r3  and  r4 
are  continuous  functions  of  0  and  8,  the  denominators  of  r3  and  r4  equal  one  when  0  and  8 
equal  zero.  The  equations  Hosking  uses  for  the  third  and  fourth  L-moments,  r3  and  r4,  in 
his  FORTRAN  search  routine  (Appendix  N)  are 


Being  a  function  of  only  two  variables,  0  and  <!>,  the  equations  for  r3  and  r4  lend 
themselves  well  to  solution  of  two  simultaneous  equations  for  two  unknowns.  Given  the  first 
four  sample  L-moments,  the  lambda  parameters  can  be  calculated  in  progression,  starting 
with  the  two  equations  above.  Once  0  and  8  have  been  determined,  the  equation  for  the 
second  L-moment  can  be  solved  directly  for  a  and  7  assuming,  as  Hosking  does,  that  0  =  7. 
Then,  the  equation  for  the  first  L-moment,  Aj,  can  be  solved  explicitly  for  £. 

4-2  Hosking’s  Solution  of  Simultaneous  Equations 

Hosking’s  algorithm  in  his  FORTRAN  subroutine,  PELTU4  (Appendix  N)  employs 
a.  form  of  the  Newton-Raphson  method  to  implicitly  solve  Equations  4.1  and  4.2  simulta¬ 
neously  for  0  and  8.  The  routine  finds  an  approximate  solution  to  the  nonlinear  system 
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F(0, 6)  =  0,  where 

F(J),6)  =&-■*)  +  (p-rO,  (4.3) 

L2 

and  an  initial  (0,S)  =  (0,0).  The  expressions  ^  and  ^  are  the  known  sample  L-moments, 
while  r3  and  r4  are  calculated  from  Equations  4.1  and  4.2 

Hosking’s  algorithm  iteratively  searches  for  the  values  of  and  6  that  yield  F(j3,6) 
within  a  pre  determined  tolerance  of  zero.  Hosking  determines  the  stepwise  changes  in  j3  and 
6  by  solving  the  derivatives  of  the  L-moment  equations,  finding  the  inverse  of  the  matrix  of 
derivatives,  and  incrementing  the  solution  vector  (/?,  6)  with  the  vector  product  of  [j^-, 
and  the  matrix  inverse. 

In  his  PELTU4  subroutine,  Hosking  chose  to  initiate  the  iterative  search  from  the 
starting  point  (/?,  <$)  =  (0,0)  because  both  f3  and  6  are  constrained  to  the  region  —1  <  x  <  1. 
Hosking  intentionally  constrains  the  parameters  to  these  values  since,  if  /?  and  8  are  unique, 
they  lie  in  this  region.  The  choice  to  begin  the  search  at  the  origin  necessitated  the  use  of 
the  revised  quantile  function  discussed  in  the  previous  section.  Hosking  allows  the  search  to 
proceed  through  twenty  iterations  before  terminating  the  algorithm  and  returning  an  error 
message  that  the  algorithm  failed  to  converge  on  an  acceptable  solution. 

4-3  The  Generalized  Lambda  Percentile  Function 

As  mentioned  earlier,  Bergevin  used  Ramberg’s  percentile  function  of  the  four- parameter 
GLD,  Equation  2.2,  to  derive  the  equations  for  Hosking’s  linear  moments  in  terms  of  the 
four  lambda  parameters.  He  derived  Equations  2.26  and  2.27  for  the  third  and  fourth  L- 
moments,  t3  and  r4  (Bergevin,  1993:  7,  63).  Like  Hosking’s  equations  for  t3  and  r4,  the 
two  equations  are  in  terms  of  only  two  variables,  A3  and  A4.  Thus,  these  equations  can  also 
be  solved  simultaneously  for  two  unknowns.  Given  the  first  four  sample  L-moments,  the 
lambda  parameters  can  be  calculated  in  progression:  Once  A3  and  A4  have  been  calculated, 
the  equation  for  the  second  L-moment  can  be  solved  directly  for  A2.  The  equation  for  the 
first  L-moment,  can  then  be  solved  explicitly  for  A4. 
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4-4  Mykytka’s  Solution  of  Simultaneous  Equations 

Mykytka  and  Ramberg  (Mykytka  and  Ramberg,  1979)  and  Mykytka  (Mykycka,  1978) 
used  non-linear  programming  methods  to  find  the  minimum  sum-of- squared  errors  between 
the  calculated  and  desired  values  of  the  third  and  fourth  conventional  moments,  o3  and  a4 

Min  /(A3,A4)  =  (or3(A3,  A4 )  —  03)“  +  (a4(A3,  A4)  —  q4)2  (4.4) 

subject  to  A3  •  A4  >  0.  (4.5) 

Expression  4.5  ensures  that  A3  and  A4  have  the  same  sign.  The  expression  in  Equation  4.4  can 
be  minimized  using  Powell’s  Algorithm  for  non-linear  function  minimization.  This  algorithm 
was  originally  designed  for  unconstrained  non-linear  function  minimization.  However,  the 
GLD  is  constrained  by  two  requirements.  First,  A3  and  A4  must  have  the  same  sign.  Second, 
the  objective  is  to  minimize  the  sum-of-squared  errors  between  the  desired  and  calculated 
values  of  o3  and  a4  subject  to  this  constraint  (see  Equations  4.4  and  4.5).  In  the  current 
implementation,  “unacceptable”  values — those  combinations  that  do  not  satisfy  Equation 
4.5 — are  eliminated  from  consideration  by  replacing  their  sum-of-squared  errors  with  a  large 
penalty.  This  penalty  is  enforced  by  defining  the  objective  function  as 

min  Z  =  /( A3,  A4) 

where 

Z  =  [(q3(A3,  A4)  —  o3)2  -f  (q4(A3,  A4)  —  a4)2];  A3  •  A4  >  0 
Z  =  10;  A3  •  A4  <  0. 

Since  an  appropriate  “match”  produces  the  desired  o3  and  a4  values  only  when  the  optimal 
value  of  the  objective  is  zero  (practically  interpreted  as  having  an  objective  function  value 
less  than  0.0002),  one  can  avoid  “unacceptable”  combinations  of  A3  and  A4  by  using  this 
penalty  (Bergevin,  1993:  18). 
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When  using  the  FORTRAN  program  given  by  Mykytka  and  Ramberg,  the  user  must 
input  an  initial  “guess”  at  the  values  for  A3  and  A4.  Since  A3  and  A4  must  have  the  same 
sign,  their  initial  values  must  be  either  both  positive  or  both  negative.  It  can  be  shown  that 
even  if  the  initial  guess  is  in  the  wrong  region  the  algorithm  will  sometimes  switch  regions 
to  find  the  optimal  solution.  For  example,  an  initial  guess  with  A3  and  A4  both  negative 
could  be  provided  when  the  optimal  solution  actually  has  them  both  positive.  But  because 
of  the  algorithm’s  search  technique,  the  possibility  of  jumping  directly  from  one  acceptable 
region  to  the  other  is  slight.  That  is,  although  the  algorithm  sometimes  does  switch  regions, 
it  usually  does  not.  If  the  user  does  not  achieve  a  satisfactory  Z-value  when  starting  in  one 
region,  he  need  simply  apply  the  algorithm  a  second  time,  changing  his  starting  point  to 
one  in  the  opposite  region  to  produce  a  better  objective  function  value  (Bergevin,  1993:  19). 
After  discovering  this,  Mykytka  revised  the  FORTRAN  algorithm  to  first  search  from  an 
initial  point  in  one  region,  (A3  =  0.05,  A4  =  0.05).  If  the  first  search  fails  to  converge  to  the 
solution,  the  algorithm  begins  another  search  from  another  region,  (A3  =  — 0.05,  A4  =  —0.05) 
to  attempt  to  reduce  the  tolerance  of  the  solution. 

A  major  part  of  Bergevin’s  thesis  was  to  study  and  test  Powell’s  Algorithm  for  its  effec¬ 
tiveness  and  robustness  when  searching  for  lambda  parameters,  especially  near  the  bound¬ 
aries  of  the  feasible  regions.  Bergevin  showed  that  Powell’s  algorithm  does  not  significantly 
limit  the  GLD.  Even  in  cases  where  the  final  A3  and  A4  values  lie  on  (or  near)  the  limiting 
boundaries,  the  algorithm  was  not  affected  by  the  penalty  function.  The  algorithm  also  per¬ 
formed  just  as  well  in  an  unconstrained  search  as  it  did  in  the  constrained  case  (Bergevin, 
’993:  32).  Mykytka’s  FORTRAN  subroutines  (Appendices  0  and  P)  compute  the  four  GLD 
parameters  from  the  first  four  conventional  moments  of  sample  data.  Because  L-moments  so 
closely  parallel  the  conventional  moments,  it  was  a  simple  matter  to  adapt  the  subroutines 
to  the  L-moment  application.  The  only  changes  needed  were  to  replace  the  conventional 
moment  equations  with  the  simpler  L-moment  equations  in  the  CALCFX  subroutine  and 
the  FN2  subroutine  (Appendix  0).  The  primary  code  for  Powell’s  algorithm  in  the  BOTM 
subroutine  (Appendix  P)  was  unchanged. 


4-5  Comparing  Newton’s  Method  and  Powell’s  Algorithm 

In  comparing  Hosking’s  Newton-Raphson  algorithm  to  Powell’s  algorithm,  three  areas 
were  examined.  Most  important  of  these  was  the  robustness  of  the  algorithm  when  searching 
for  lambda  parameters  near  the  limits  of  the  feasible  regions.  Second  was  the  accuracy  of 
the  lambda  parameters  computed  for  the  sample  L-moments.  The  final  area  was  the  speed 
and  efficiency  of  the  computer  search  subroutines. 

4-5.1  Robustness.  The  most  important  criteria  for  a  search  method  is  its  ability 
to  converge  to  a  solution.  A  search  method  will  most  likely  encounter  difficulty  when  it 
approaches  the  boundary  of  the  feasible  region  or  a  point  of  discontinuity  of  a  function. 
Both  Newton’s  method  and  Powell’s  algorithm  must  deal  with  boundary  conditions  as  the 
parameters  approach  1  or  —1.  Bergevin’s  equations  also  contain  a  discontinuity  when  the 
lambda  parameters  approach  zero. 

While  the  Newton-Raphson  algorithm  is  a  “classic”  numerical  method,  Powell’s  al¬ 
gorithm  benefits  from  Mykytka’s  and  Bergevin’s  extensive  verification  and  validation  ef¬ 
forts  for  this  particular  application.  Powell’s  algorithm  is  better  for  dealing  with  discon¬ 
tinuous  and  multi-modal  objective  functions  (Bergevin,  1993:  29),  and  is  appropriate  to 
Bergevin’s  discontinuous  L-moment  equations.  In  addition,  during  comparison  tests  of  the 
two  methods,  the  Newton-Raphson  method  sometimes  converged  on  an  unacceptable  solu¬ 
tion  (A3  >  0,  A4  <  0)  for  Gamma  distributions.  Newton’s  method  also  failed  to  converge 
on  a  solution  for  some  Beta  and  Exponential  distributions.  Both  of  these  distributions  have 
parameter  values  which  lie  near  the  region  boundaries.  For  these  reasons,  Powell’s  algorithm 
is  considered  to  be  more  robust  than  the  Newton-Raphson  method. 

4-5.2  Accuracy.  To  test  the  accuracy  of  each  algorithm,  a  FORTRAN  subroutine 
was  developed  to  re-compute  the  L-moments,  given  the  estimated  lambda  parameters  (Ap¬ 
pendix  Q).  When  comparing  the  re-computed  L-moments  to  the  original  sample  L-moments, 
Powell's  algorithm  was  consistently  accurate  to  six  decimal  places  on  all  four  L-moments. 


Newton’s  method  was  also  accurate  to  six  decimal  places  when  computing  the  first  and  sec¬ 
ond  L-moments.  But,  it  was  accurate  to  only  four  decimal  places,  and  occasionally  accurate 
to  only  two  decimal  places  on  the  third  and  fourth  L-moments.  But  when  the  lambda  pa¬ 
rameters  from  both  routines  were  used  to  plot  the  PDF’s  of  their  respective  distributions, 
there  was  no  visible  difference  between  the  plots.  Therefore,  the  results  provided  by  the 
Newton-Raphson  method  and  Powell’s  algorithm  were  equally  accurate. 

4-5.3  Speed  and  Efficiency.  The  Newton-Raphson  method  benefits  from  simplicity 
and  quadratic  convergence.  In  tests,  Hosking’s  routine  typically  converged  to  an  answer  in 
less  than  ten  iterations,  for  all  types  of  distributions.  This  fast  rate  of  convergence,  combined 
with  a  small  number  of  instructions,  makes  the  Newton-Raphson  method  quite  fast  and 
efficient.  Powell’s  algorithm  was  technically  less  efficient,  sometimes  requiring  sixty  iterations 
to  converge  to  a  solution.  This,  coupled  with  about  three  times  as  many  instructions  per 
iteration  as  the  Newton-Raphson  routine,  makes  Powell’s  algorithm  slower  and  less  efficient. 
Still,  the  overall  run  time  for  either  routine  is  a  matter  of  microseconds  of  CPU  time. 

Both  routines  are  part  of  a  program  which  runs  in  real-time  and  interacts  with  a  human 
operator  arid  batch  files  stored  in  memory.  In  this  type  of  application,  both  input /output 
to  memory  and  wait  times  from  screen  prompt  to  keyboard  input  are  long  compared  to 
actual  CPU  time  spent  executing  instructions.  During  testing,  there  was  no  discernible 
delay  between  the  last  keyboard  input  and  the  screen  prompt  for  the  next  input,  using 
either  algorithm.  So,  while  the  Newton-Raphson  method  is  technically  faster  than  Powell’s 
algorithm,  they  are  practically  both  instantaneous.  Therefore,  there  is  no  discrimination 
between  the  two  based  on  speed  and  efficiency. 

In  summary,  Powell’s  algorithm  is  more  robust  than  the  Newton-Raphson  method. 
Both  methods  are  equally  suitable  with  respect  to  accuracy,  speed,  and  efficiency.  But,  since 
robustness  is  the  most  important  quality,  Powell’s  algorithm  is  the  overall  best  choice  for  a 
computer  search  routine. 


V.  Design  of  the  Monte  Carlo  Experiment 


5. 1  Monte  Carlo  Experiments 

A  Monte  Carlo  simulation  study  was  done  to  compare  the  effects  of  sampling  variability 
in  the  method  of  linear  moments,  the  conventional  method  of  moments,  and  the  alternative 
method  of  moments  using  Q3  and  Q4  as  measures  of  symmetry  and  tailweight.  Henceforth, 
these  methods  are  referred  to  as  linear  moments,  conventional  moments,  and  alternate  mo¬ 
ments,  respectively.  Each  of  the  three  methods  was  used  to  estimate  the  lambda  parameters 
from  nine  different  sets  of  thirty  random  samples  each.  The  nine  sets  were  combinations  of 
three  different  Lambda  distributions  approximating  the  Normal,  Gamma,  and  Exponential 
distributions,  respectively;  and  three  different  sample  sizes:  n  =  25,  n  =  50,  and  n  =  100. 
Random  samples  were  generated  using  the  GLD  quantile  function  to  transform  uniform  [0,1] 
pseudo-random  variates  into  variates  from  the  selected  distribution.  The  pseudo-random 
number  generator  DURAND,  developed  by  IBM,  was  employed  and  is  described  in  Appendix 
K.  The  same  thirty  seed  values  were  used  to  generate  the  thirty  samples  for  each  of  the  nine 
sample  sets  (the  seed  values  used  were  1010,  1020,  1030,  ...  ,  1300). 

The  sampling  variability  effects  were  compared  by  plotting  -  for  each  method,  distri¬ 
bution,  and  sample  size  -  the  estimated  probability  density  functions  (PDF)  for  the  first 
twenty-five  of  the  thirty  replications  on  the  same  graph.  The  software  package  used  for 
plotting  the  PDF’s  limited  the  number  of  simultaneous  plots  to  twenty-five.  Kolmogorov- 
Smirnov  goodness-of-fit  statistics  were  also  collected  from  each  replication  by  comparing 
the  fitted  cumulative  distribution  functions  (CDF)  to  that  of  the  underlying  theoretical 
distribution  from  which  the  sample  was  randomly  generated. 

The  overall  goal  of  this  study  was  to  determine  whether  or  not  the  method  of  linear 
moments  is  less  susceptible  to  sampling  variability  effects  than  the  method  of  conventional 
moments.  Therefore,  the  results  of  the  Monte  Carlo  experiments  using  conventional  moments 
were  used  as  a  “control  group,”  and  the  results  of  the  other  experiments  were  compared  to 


these  baseline  experiments.  Thus,  the  conclusions  drawn  were  expressed  in  relative,  rather 
than  absolute  measures. 

5.2  Data  Sample  Sizes 

Because  the  goal  is  to  improve  the  GLD’s  ability  to  estimate  the  parameters  of  the 
underlying  distributions  from  small  samples  and  reduce  the  sampling  variability  effect,  it 
is  necessary  to  experiment  with  “small”  sample  sizes.  Mykytka  (1978:60)  had  previously 
determined  that  conventional  moments  were  adequate  for  sample  sizes  as  small  as  150,  and 
that  smaller  sample  sizes  introduced  unacceptable  sampling  variability.  Therefore,  sample 
sizes  of  n  =  25,  n  =  50,  and  n  =  100  were  selected  as  representative  “small”  sample  sizes. 
These  sample  sizes  are  also  more  in  line  with  samples  being  used  in  today’s  modeling  and 
testing  communities  because  of  cost  and  time  constraints. 

5.3  Distributions  Selected 

In  order  to  get  a  representative  cross-section  of  various  distribution  types  and  shapes, 
three  Lambda  distributions  whose  PDF’s  approximate  the  shapes  of  Normal,  Gamma,  and 
Exponential  PDF’s,  respectively,  were  chosen.  These  three  distributions  are  commonly-used 
and  well-known  to  anyone  familiar  with  statistics.  As  seen  in  Figure  2.2,  thev  represent 
a  broad  range  of  possible  combinations  of  skewness  and  kurtosis.  These  three  Lambda 
distributions  were  chosen  to  have  the  same  first  four  moments  as  their  respective  Normal, 
Gamma,  and  Exponential  counterparts.  To  remind  the  reader  that  they  are  the  Lambda 
approximations  to  these  distributions,  they  will  be  referred  to  with  their  names  in  quotes, 
i.e.,  as  the  “Normal,”  “Gamma,”  and  “Exponential”  distributions,  respectively.  Table  5.1 
summarizes  the  characteristics  of  the  chosen  distributions. 


Parameter 

Approximate  Distribution  Shape 

“Normal” 

“Gamma” 

“Exponential” 

Mean 

0 

0.8 

1 

Variance 

1 

1 

1 

Skewness 

0 

v/2 

2 

Kurtosis 

3 

6 

9 

Ai 

0 

0 

0 

A2 

0.1975 

0.04134 

-0.001632 

A3 

0.1349 

0.005674 

-0.9159  x  10-5 

A4 

0.1349 

0.04046 

-0.001621 

Ai 

0 

7.6558  x  10~5 

6.0309  x  10“5 

a2 

0.56382 

0.52905 

0.50065 

r3 

0 

0.23425 

0.33031 

74 

0.12447 

0.15455 

0.16723 

M 

0 

0.8 

1 

a2 

1 

1 

1 

Qs 

1 

2.64S1 

4.4743 

Q< 

2.5959 

2.7856 

2.8676 

Table  5.1  Moments  and  parameters  of  distributions  used  in  the  experiment 


5.4  Methods  of  Fit 

The  purpose  of  this  study  was  to  determine  whether  or  not  the  method  of  linear 
moments  is  a.  better  method  for  fitting  a  distribution  to  sample  data  than  the  method  of 
conventional  moments.  Obviously,  these  two  methods  were  included  in  the  experiment.  The 
alternate  method  of  fit,  using  Q 3  and  Q4,  was  also  be  included  in  the  experiment.  As  Mykytka 
(1978)  has  shown,  the  method  of  alternate  moments  suffers  less  from  sampling  variability 


effects  than  conventional  moments.  In  addition,  alternate  moments  use  a  similar  method 
for  computing  lambda  parameters.  Like  conventional  moments  and  linear  moments,  it  uses 
Powell’s  search  algorithm  to  find  the  third  and  fourth  lambda  parameters.  The  objective 
function  subrout  ine  is  the  only  difference  between  the  three  methods. 

Although  Chou  (1988)  showed  that  the  method  of  maximum  likelihood  estimators 
(MLE’s)  produces  better  fits  and  better  parameter  estimates  than  conventional  moments, 
it  was  not  included  in  this  experiment  since  its  cpu-time  requirements  remain  excessive  for 
most  applications. 

5.5  Measures  of  Effectiveness 

5.5. 1  Statistical  Summary.  Table  5.2  is  an  example  of  the  summary  of  information 
collected  from  each  experiment.  Each  experiment  consisted  of  fitting  the  GLD  to  thirty 
different  data  samples  using  each  of  the  three  estimation  methods.  The  notes  at  the  bottom 
of  the  table  describe  the  experiment.  The  fourteen  summary  statistics  listed  along  the  left- 
hand  side  of  the  table  were  collected  from  each  replication  of  the  experiment.  For  each  of 
the  fourteen  categories,  the  minimum  and  maximum  values  of  each  statistic  over  the  thirty 
replications  are  displayed  in  the  table.  In  addition,  the  mean  and  variance  were  calculated  for 
the  group  of  thirty  values  in  each  category  and  were  also  displayed  in  the  table.  The  rows 
labelled  “Theoretical  K-S  Statistics,”  and  “Empirical  K-S  Statistics”  contain  information 
about  the  Kolmogorov-Smirnov  statistics  collected  from  comparisons  of  the  theoretical  CDF 
to  the  fitted  CDF,  and  of  the  empirical  CDF  to  the  fitted  CDF,  respectively. 


Statistic 


Experiment 


Result 

Minimum 

Maximum 

Average 

Variance 

1st 

-0.32920000 

0.33291600 

-0.08192267 

0.03075791 

Moments  of 

2nd 

0.46495600 

0.00435006 

Sample  Data 

3rd 

-0.14980300 

0.16846000 

0.02440728 

0.00632024 

4th 

-0.01568870 

0.00315052 

'V 

-2.07855000 

1.39302000 

-0.28595313 

0.73005042 

Lambda 

A  2 

0.22071600 

0.49963800 

0.19608461 

0.02915614 

Parameters 

A3 

-0.08825540 

1.00248000 

0.18008392 

0.06060880 

a4 

-0.12274100 

1.11750000 

0.23446677 

0.07488433 

Theoretical 

MIN 

0.00000211 

0.00132966 

0.00015583 

0.00000008 

K-S 

AVG 

0.01000620 

0.05861190 

0.03572510 

0.00016930 

Statistics 

MAX 

0.02489270 

0.18884700 

0.10145870 

0.00149412 

Empirical 

MIN 

0.00000000 

0.00024583 

0.00002845 

0.00000000 

K-S 

AVG 

0.01290350 

0.02903130 

0.02055428 

0.00002040 

Statistics 

MAX 

0.06148730 

0.14373100 

0.09572504 

0.00051787 

Experiment  Number:  1A 


Samples:  30  samples  of  25  elements  each 

Original  function:  Lambda  approximation  of  a  Normal  function 
Method  of  fit:  method  of  linear  moments 


Table  5.2  Statistical  summary  of  Monte  Carlo  experiment  results 


5.5.2  PDF’s  and  Visual  Fits.  Figures  5.1,  5.2,  and  5.3  are  examples  of  the  three 
distributions  used  in  the  experiment.  Each  figure  contains  both  a  plot  of  the  theoretical 
density  function  and  a  plot  of  the  density  function  estimated  from  a  1000-element  sample  by 
the  method  of  linear  moments.  The  experimental  results  in  Appendices  B  through  .J  contain 
plots  of  25  estimated  density  functions  from  their  respective  sets  of  data  samples.  A  com  par- 


Figure  5.1  Comparison  of  the  theoretical  “Normal”  PDF  to  the  PDF  fit  from  a  sample  of 
1000  observations 

ison  of  the  different  plots  shows  the  methods’  relative  abilities  to  consistently  estimate  the 
parameters  of  the  underlying  distributions.  These  plots  are  suitable  for  visual  comparisons 
only,  and  will  be  discussed  in  Chapter  Six.  No  effort  is  made  to  collect  quantitative  data 
from  these  graphs. 

5.5.3  CDF’s  and  Kolmogorov-Smirnov  Statistics.  Figure  5.4  is  an  example  of  the 
plot  of  the  CDF  of  the  theoretical  “Normal”  distribution  and  a  plot  of  the  CDF  of  a  distri¬ 
bution  estimated  from  a  25-element  sample  by  the  method  of  linear  moments.  Figure  5.5  is 
an  example  of  the  plot  of  the  empirical  CDF  of  a  25-element  sample  and  the  CDF  plot  of  the 
distribution  estimated  from  that  sample  by  the  method  of  L-moments.  From  such  informa¬ 
tion,  various  Kolmogorov-Smirnov  goodness-of-fit  statistics  were  computed.  In  particular, 
using  the  “Normal”  distribution  shown  in  Figure  5.4  as  an  example,  the  vertical  distance 
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0.6, 


Figure  5.2  Comparison  of  the  theoretical  “Gamma”  PDF  to  the  PDF  fit  from  a  sample  of 
1000  observations 

between  the  two  CDFs  were  measured  at  598  evenly-spaced  points  along  the  horizontal  axis, 
beginning  at  -2.99  and  ending  at  2.99.  From  these  measurements,  the  minimum  distance, 
maximum  distance,  and  the  average  distance  were  recorded.  These  distances  are  referred 
to  as  the  minimum,  maximum,  and  average  K-S  values.  A  comparison  of  a  fitted  CDF  to 
a  theoretical  CDF  attempts  to  capture  the  method’s  ability  to  estimate  the  underlying  dis¬ 
tribution  of  the  population.  The  maximum  Kolmogorov-Smirnov  value  typically  represents 
a  worst-case  error  between  the  fitted  and  theoretical  CDF’s.  The  average  Kolmogorov- 
Smirnov  value  represents  a  measure  of  overall  success  in  estimating  a  fit  to  the  respective 
population  or  sample. 

The  results  of  applying  these  measures  of  fit.  to  the  results  of  the  Monte  Carlo  experi¬ 
ment  are  described  in  Chapter  Six. 
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Figure  5.3  Comparison  of  the  theoretical  “Exponential”  PDF  to  the  PDF  fit  from  a  sample 
of  1000  observations 
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Fitted  CDF  vs.  Theoretical  CDF 


ure  5.4  Comparison  of  the  theoretical  “Normal”  CDF  to  the  CDF  fit  from  a  sample 
25  observations 


Fitted  CDF  vs.  Empirical  CDF 


Figure  5.5  Comparison  of  the  empirical  CDF  of,  and  a  CDF  fit  to,  a  sample  of  25  observa¬ 
tions  from  the  “Normal”  distribution 
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VI.  Results  of  the  Monte  Carlo  Experiments 

6. 1  Tables  of  Experimental  Results 

Tables  6.1  through  6.10  and  Figures  6.1  through  6.9  summarize  the  data  collected 
from  the  Monte  Carlo  simulation  study.  The  data  is  grouped  to  compare  the  three  different 
methods  of  fit.  The  data  from  the  experiments  can  be  found  in  Appendices  B  through  J. 

6.2  Fitted  Distributions  vs.  the  Underlying  Theoretical  Distribution 

6.2.1  Comparing  Probability  Density  Function  Plots.  Figures  6.1  through  6.9 
display  the  underlying  theoretical  probability  density  functions  of  the  sample  data  used  in 
the  Monte  Carlo  experiment.  Alongside  the  theoretical  PDF  are  twenty-five  examples  of 
PDF’s  fitted  to  sample  data  by  each  of  the  three  methods  of  moments.  Appearing  clockwise 
from  the  theoretical  PDF  are  the  linear,  conventional,  and  alternate  methods,  respectively. 
The  same  data  samples  are  used  by  each  of  the  three  estimation  methods. 

The  comparison  in  Figure  6.1  of  PDF’s  fitted  from  25-element  samples  shows  that  a 
substantial  amount  of  ‘‘sampling  variability,”  or  “noise,”  is  transferred  from  the  sample  data 
to  the  fitted  distributions  by  all  three  methods.  Note  that  each  set  of  examples  contains  two 
remarkably  poor  fits.  An  examination  of  the  individual  samples  shows  that  the  “poor  fits” 
are  caused  by  the  same  two  samples  in  each  case.  Further,  these  two  fits  have  the  highest 
Kolmogorov-Smirnov  statistics  found  in  the  entire  set  of  thirty  samples. 

The  comparison  in  Figure  6.2  of  PDF’s  fitted  from  50-element  samples  shows  the  noise 
carried  over  from  the  sample  data  is  reduced  by  the  larger  sample  size,  but  suggests  that 
the  different  methods  of  fit  remain  basically  indistinguishable.  Again,  note  that  each  set 
contains  one  poorly-fitted  distribution.  Each  is  a  result  of  the  same  data  sample.  That 
sample  also  produces  the  highest  K-S  value  in  the  set  of  samples. 

The  comparison  in  Figure  6.3  of  PDF’s  fitted  from  100-element  samples  shows  that  the 
noise  is  substantially  reduced  by  the  sample  size  for  all  three  methods  of  moments  and  that 
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Figure  6.1  Theoretical  “Normal”  PDF  vs.  Methods  of  Fit  of  25-Element  Samples 

there  is  no  visible  difference  between  groups  of  distributions  fitted  by  the  different  methods. 
Also,  this  sample  size  was  apparently  sufficient  to  eliminate  any  radically-differcnt  fits  caused 
by  sampling  variability. 

The  comparison  in  Figure  6.4  of  PDF's  fitted  from  25-element  samples  selected  from 
the  "Gamma"  distribution  shows  that  a  significant  amount  of  noise  is  again  carried  to  the 
estimated  PDF  s  by  all  three  methods.  Note  that,  the  conventional  method  of  moments 
lias  noticeable  more  'poor  fits  than  t  he  other  two  methods.  Several  fits  have  taken  on  an 
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Figure  6.2  Theoretical  “Normal”  PDF  vs.  Methods  of  Fit  of  50-Element  Samples 


Exponential-like  “ski  slope”  shape,  rather  than  the  desired  skewed  unimodal  shape  of  the 
“Gamma”  distribution. 

The  comparison  in  Figure  6.5  of  PDF’s  fitted  from  50-element  samples  shows  that  the 
three  methods  again  exhibit  roughly  equal  levels  of  sampling  variability,  and  that  the  effects 
of  sampling  variability  are  reduced  by  the  increase  in  sample  size. 

The  comparison  in  Figure  6.6  of  PDF's  fitted  from  100-element  samples  shows  that 
the  increase  in  sample  size  continues  to  reduce  the  effects  of  sampling  variability  for  all  three 


Figure  6.3  Theoretical  “Normal”  PDF  vs.  Methods  of  Fit  of  100-Element  Samples 

methods  of  fit.  Note,  however,  that  this  sample  size  no  longer  eliminates  “poor  fits”  by 
conventional  moments. 

The  comparison  in  Figure  6.7  of  PDF’s  fitted  from  25-element  samples  shows  that,  for 
the  “Exponential”  case,  the  different  methods  of  fit  show  different  levels  of  resistance  to  the 
effects  of  sampling  variability.  The  conventional  moments  show  the  least  resistance,  with 
the  largest  group  of  “Exponential  shaped  distributions,  while  the  linear  moments  show  the 
most  resistance,  with  the  “tightest”  group  of  “Exponential”  shaped  distributions. 
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Figure  6.4  Theoretical  “Gamma”  PDF  vs.  Methods  of  Fit  of  25-Element  Samples 

The  comparison  in  Figure  6.8  of  PDF’s  fitted  from  50-element  samples  shows  that  the 
sample  size  has  reduced  the  effective  noise  of  the  sample,  but  the  relative  “behavior”  of  the 
different  methods  shown  in  the  previous  Figure  still  remains  in  effect. 

The  comparison  in  Figure  6.9  of  PDF’s  fitted  from  100-element  samples  shows  that 
this  sample  size  is  insufficient  in  this  case  to  mask  the  noise  in  sample  data  from  any  of  the 
methods  of  fit.  This  suggests  that  as  the  skewness  and  kurtosis  of  an  underlying  distribution 
increases,  the  size  of  sample  data  must  also  increase  to  mask  a  given  level  of  sampling 
variability.  The  relative  abilities  of  the  different  methods  to  reduce  the  sampling  variability 
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Figure  6.5  Theoretical  “Gamma”  PDF  vs.  Methods  of  Fit  of  50-Element  Samples 

effects  remains  visible.  The  linear  moments  remain  most  effective,  while  the  conventional 
moments  are  still  least  effective. 

This  visual  comparison  of  groups  of  fitted  distributions  does  provide  a  subjective  feel 
for  which  methods  of  fit  are  better  or  worse,  but  it  does  not  provide  an  objective  or  more 
rigorous  comparison.  In  order  to  make  a  quantitative  assessment  of  the  relative  differences 
in  the  methods  of  fit,  an  effort  was  made  to  determine  the  percentage  of  fitted  distributions 
which  accurately  capture  the  shape  of  the  underlying  distribution. 
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Figure  6.6  Theoretical  “Gamma”  PDF  vs.  Methods  of  Fit  of  100-Element  Samples 

6.2.2  Summarizing  Probability  Density  Function  Plots.  Table  6.1  summarizes  the 
observations  made  from  the  visual  comparisons  of  the  PDF’s  estimated  from  the  Monte 
Carlo  simulation  study  and  lists  the  number  of  “poor  fits”  observed  in  each  set  of  25  PDF’s 
estimated  by  the  linear,  conventional,  and  alternate  methods,  respectively.  Any  plotted 
curve  that  does  not  resemble  the  shape  of  the  underlying  distribution  constitutes  a  poor 
fit.  Figures  6.10  through  6.12  provide  examples  of  fitted  PDF’s  which  do  not  resemble  the 
shapes  of  the  theoretical  “Normal,”  “Gamma,”  and  “Exponential”  PDF’s,  respectively. 
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Figure  6.7  Theoretical  “Exponential”  PDF  vs.  Methods  of  Fit  of  25-Element  Samples 

A  distribution  fitted  from  “Normal”  distributed  sample  data  would  be  considered  a 
“poor  fit”  if  it  was  not  unimodal,  roughly  symmetrical,  and  approximately  “bell-shaped.” 
A  distribution  fitted  from  “Gamma”  distributed  sample  data  would  be  considered  a  “poor 
fit”  if  it  was  not  unimodal,  asymmetrical,  and  skewed  to  the  left.  A  distribution  fitted  from 
“Exponential”  distributed  sample  data  would  be  considered  a  “poor  fit”  if  it  was  not  concave, 
asymmetrical  and  skewed  to  the  left,  and  approximately  “ski  slope-shaped.” 

Table  6.1  shows  that  all  three  methods  of  fit  are  very  capable  of  estimating  good  fits  to 
the  “Normal”  distribution,  at  any  of  the  tested  sample  sizes.  The  methods  of  linear  moments 
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Figure  6.8  Theoretical  “Exponential”  PDF  vs.  Methods  of  Fit  of  50-Element  Samples 

and  alternate  moments  were  noticeably  better  than  conventional  moments  at  estimating  the 
PDF’s  of  samples  from  the  “Gamma”  distribution.  While  the  “Exponential”  distribution 
was  the  most  difficult  to  estimate  for  all  three  methods,  it  was  particularly  difficult  for  the 
method  of  conventional  moments.  That  rriethod  was  effectively  unreliable  for  estimating 
a  good  fit  from  an  “Exponential”  distributed  sample.  The  method  of  linear  moments  was 
better  than  the  method  of  alternate  moments,  but  not  to  the  degree  of  its  advantage  over 
the  method  of  conventional  moments.  Also,  for  a  given  sample  size,  sample  data  from 
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Figure  6.9  Theoretical  “Exponential”  PDF  vs.  Methods  of  Fit  of  100-Element  Samples 

distributions  with  higher  skewness  and  kurtosis  values  transferred  more  noise  to  the  fitted 
distributions  through  all  three  methods  of  moments. 

Although  this  summary  of  the  fitted  PDF’s  provides  an  intuitive  feel  for  the  relative 
advantages  or  disadvantages  of  the  various  methods  of  moments,  it  is  still  rather  subjective. 
In  addition,  while  the  assessments  made  so  far  are  useful  in  comparing  different  methods  to 
each  other,  they  are  not  so  useful  for  comparing  them  against  the  theoretical  distribution 
itself. 
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Figure  6.10  Examples  of  “poor  fits”  to  the  theoretical  “Normal”  PDF 


The  Kolmogorov-Smirnov  goodness-of-fit  test  is  a  quantitative  measure  of  the  accu¬ 
racy  with  which  a  fitted  CDF  matches  the  underlying  CDF.  Various  K-S  statistics  were 
collected  from  the  fitted  CDF  of  every  data  sample  and  compiled  to  provide  a  quantitative 
measure  of  the  goodness-of-fit  of  each  of  the  three  methods  of  moments,  both  for  comparison 
to  the  theoretical  distribution,  and  to  each  other. 

6.2.3  Summarizing  the  Kolmogorov-Smirnov  Statistics.  Table  6.2  summari^  the 
Kolmogorov-Smirnov  (K-S)  statistics  collected  from  each  of  the  27  Monte  Carlo  experiments. 
The  “Average  K-S”  Column  lists  the  average  over  the  30  samples  of  the  average  difference 
between  the  estimated  and  theoretical  distribution  functions.  The  “Maximum  K-S”  Column 
lists  the  average  over  the  30  samples  of  the  maximum  difference  between  the  same  functions. 
The  “Max  K-S  Variance”  Column  lists  the  variance  of  the  30  values  which  form  the  average 
in  previous  column. 
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Figure  6.11  Examples  of  “poor  fits”  to  the  theoretical  “Gamma”  PDF 


When  comparing  the  different  methods  of  moments  for  like  sample  distributions  and 
sample  sizes,  there  is  little  difference  among  the  averages  of  the  K-S  statistics.  The  K-S 
averages  of  the  linear  and  alternate  methods  are  identical  to  one  significant  digit,  and  in  most 
cases,  identical  to  two  significant  digits.  The  conventional  method  is  only  slightly  worse  for 
the  “Gamma”  and  “Exponential”  distributions.  Based  on  these  averages,  no  method  shows 
an  obvious  advantage  over  another. 

However,  there  are  noticeable  differences  in  the  variances  of  these  sample  sets  among 
the  different  methods  of  fit.  For  the  “Normal”  distribution,  all  three  methods  perform  equally 
well  for  all  sample  sizes.  But  for  the  “Gamma”  and  “Exponential”  distributions,  particularly 
for  the  smaller  sample  sizes,  the  linear  moments  and  alternate  moments  achieve  smaller 
variances  than  the  conventional  moments.  Both  methods  average  a  noticeable  decrease 
in  variance  over  the  conventional  moments  for  all  sizes  of  “Gamma”  and  “Exponential” 


6-12 


Figure  6.12  Examples  of  “poor  fits”  to  the  theoretical  “Exponential”  PDF 

distributed  samples.  This  reduced  variance  shows  that  the  methods  of  linear  and  alternate 
moments  produce  more  consistent  estimates  for  skewed  distributions. 
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Number  of  “Poor  Fits”  out  of  ‘25  Samples 

Distribution 

Sample 

Linear 

Conventional 

Alternate 

Shape 

Size 

Moments 

Moments 

Moments 

25 

2 

3 

3 

“Normal” 

50 

1 

1 

1 

100 

0 

1 

0 

25 

2 

14 

4 

“Gamma” 

50 

3 

9 

3 

100 

0 

7 

0 

25 

8 

12 

7 

“Exponential” 

50 

9 

15 

6 

100 

4 

14 

3 

Table  6.1  Summary  of  PDF  Plots 


Distribution 

Shape 

Sample 

Size 

Method 

of  Moments 

Average 

K-S 

Maximum 

K-S 

Max  K-S 

Variance 

linear 

0.03573 

0.10150 

.00149 

25 

conventional 

0.03548 

0.09837 

.00104 

alternate 

0.03524 

0.09806 

.00119 

linear 

0.02452 

0.06728 

.00086 

“Normal” 

50 

conventional 

0.02522 

0.06705 

.00100 

alternate 

0.02454 

0.06675 

.00087 

1 

inear 

0.01886 

0.05100 

.00051 

100 

conventional 

0.01836 

0.04879 

.00048 

alternate 

0.01863 

0.05116 

.00056 

linear 

0.02790 

0.09577 

.00132 

25 

conventional 

0.02837 

0.09355 

.00104 

alternate 

0.02824 

0.09368 

.00099 

linear 

0.02040 

0.06530 

.00077 

“Gamma” 

50 

conventional 

0.02319 

0.07323 

.00092 

alternate 

0.02088 

0.06534 

.00081 

linear 

0.01549 

0.05020 

.00044 

100 

conventional 

0.01741 

0.05959 

.00090 

alternate 

0.01555 

0.04911 

.00045 

1 

inear 

0.02480 

0.09237 

.00128 

25 

conventional 

0.02693 

0.11733 

.00156 

alternate 

0.02529 

0.08874 

.00174 

1 

inear 

0.01895 

0.06403 

.00063 

“Exponential” 

50 

conventional 

0.02239 

0.09521 

.00112 

alternate 

0.01868 

0.05966 

.00082 

inear 

0.01450 

0.04547 

.00038 

100 

conventional 

0.01663 

0.07962 

.00096 

alternate 

0.01428 

0.04437 

.00043 

Table  6.2  K-S  Statistic  Summary  of  fitted  CDF’s  vs.  the  theoretical  CDF 
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Average  K-S  Values 
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Figure  6.13  Average  K-S  Values  of  fitted  CDF’s  vs.  the  theoretical  CDF 

Figures  6.13  through  6.15  display  the  information  from  the  “Average  K-S,”  “Maximum 
K-S,”  and  “Max  K-S  Variance”  columns  of  Table  6.2,  respectively.  These  are  helpful  in 
highlighting  the  differences  between  the  methods. 

Figure  6.13  clearly  shows  that  the  linear  and  alternate  moments  are  comparable  while 
the  conventional  moments  produce  slightly  higher  averages  of  the  K-S  values.  There  are  two 
downward  trends  depicted  in  the  chart.  The  average  K-S  values  decrease  as  the  sample  size 
increases,  and  as  the  skewness  and  kurtosis  of  the  theoretical  distribution  increase. 

Figure  6.14  shows  even  more  disparity  between  the  performance  of  the  method  of 
conventional  moments  and  the  other  two  methods.  Note  that  for  all  cases  of  “Normal  ”  dis¬ 
tributed  sample  data,  all  three  methods  perform  equally.  For  “Gamma”  and  “Exponential” 
distributed  sample  data,  the  K-S  values  associated  with  the  conventional  method  of  moments 
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Maximum  K-S  Values 
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Figure  6.14  Maximum  K-S  Values  of  fitted  CDF’s  vs.  the  theoretical  CDF 

become  increasingly  larger  than  those  for  the  other  moments  as  the  skewness  and  kurtosis 
increases. 

Figure  6.15  suggests  some  instability  in  the  methods  at  low  sample  sizes  as  indicated 
by  by  the  unusual  behavior  shown  in  the  chart  for  25-element  sample  data.  For  sample 
sizes  of  fifty  elements  and  up,  the  previous  trends  continue.  That  is.  the  linear  and  alternate 
methods  perform  about  equally,  while  the  conventional  moments  do  not  perform  as  well.  The 
variance  in  K-S  values  also  tends  to  decrease  as  the  sample  size  and  skewness  and  kurtosis 
of  the  distribution  increase. 

Table  6.2  and  Figures  6.13  through  6.15  provide  a  quick  overall  picture  of  the  objective 
measures  of  the  goodness-of-fit  for  each  of  the  three  methods  of  fit.  The  mean  and  variance 
of  the  K-S  statistics  provide  an  idea  of  which  method  will  perform  well  in  the  long  run. 
However,  information  on  which  method  will  most  often  provide  the  closest  fit  is  also  useful. 
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Variances  of  Maximum  K-S  Values 
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Figure  6.15  Variances  of  Maximum  K-S  Values  of  fitted  CDF’s  vs.  the  theoretical  CDF 

6.2.4  Kolmogorov-Smirnov  Statistics  of  Individual  Data  Samples.  Tables  6.3  and 
6.4  list  the  number  of  times  that  each  method  of  fit  scored  the  lowest  K-S  value  among  each 
in  a  set  of  thirty  samples.  The  tables  are  divided  into  summaries  of  maximum  and  average 
K-S  values.  Each  table  is  further  broken  down  by  the  underlying  distribution  and  size  of  the 
samples  in  the  set. 
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Fitted  Distribution 

vs.  Theoretical  Distribution 

MAX  K-S  of 

MAX  K-S  of 

MAX  K-S  of 

Distribution  Shape 

Sample  Size 

Linear 

Method 

Conventional 

Method 

Alternate 

Method 

25 

9 

12 

9 

“Normal” 

50 

S 

12 

10 

100 

14 

12 

4 

25 

13 

7 

10 

“Gamma” 

50 

13 

8 

9 

100 

13 

6 

11 

25 

11 

7 

12 

“Exponential” 

50 

11 

4 

15 

100 

16 

3 

11 

Table  6.3  Number  of  Lowest  MAX  K-S  Theoretical  Differences  per  set  of  30  samples 


Table  6.3  shows  that  the  method  of  conventional  moments  most  often  provides  the 
best  fit  to  the  underlying  theoretical  “Normal”  distribution.  The  method  of  linear  moments 
most  often  produces  the  closest  fit  to  the  underlying  theoretical  “Gamma”  distribution. 
Both  the  linear  and  alternate  moments  are  good  for  estimating  the  underlying  theoretical 
“Exponential”  distribution.  None  of  the  methods  show  a  distinct  superiority  as  the  sample 
size  is  reduced.  The  lack  of  a  consistently  superior  method  suggests  a  lack  of  correlation 
between  the  method  of  fit  used  on  sample  data  and  the  quality  of  fit  in  general. 
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Fitted  Distribution 

vs.  Theoretical  Distribution 

AVG  K-S  of 

AVG  K-S  of 

AVG  K-S  of 

Distribution  Shape 

Sample  Size 

Linear 

Method 

Conventional 

Method 

Alternate 

Method 

25 

10 

8 

12 

“Normal” 

50 

8 

13 

9 

100 

14 

10 

6 

25 

15 

7 

8 

“Gamma” 

50 

17 

6 

7 

100 

14 

7 

9 

25 

12 

12 

6 

“Exponential” 

50 

11 

9 

10 

100 

12 

6 

12 

Table  6.4  Number  of  Lowest  AVG  K-S  Theoretical  Differences  per  set  of  30  samples 


Table  6.4  shows  that  the  methods  of  linear  and  conventional  moments  estimate  the 
theoretical  “Normal”  distribution  well.  The  method  of  linear  moments  is  best  at  fitting  both 
the  theoretical  “Gamma”  and  “Exponential”  distribution  shapes.  Again,  no  method  shows 
a  clear  superiority  for  smaller  sample  sizes  across  all  distributions  (although  L-moments  look 
at  least  as  good  as  the  others). 

Having  examined  the  distributions  visually  and  having  looked  at  quantitative  goodness- 
of— fit  data  taken  from  the  distributions,  another  source  of  information  is  the  set  of  lambda 
parameters  estimated  from  the  sample  data.  Comparing  the  parameter  estimates  to  their 
theoretical  values  provides  another  measure  of  a  method’s  ability  to  determine  a  distribution 
from  sample  data. 
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6.2.5  Summarizing  Lambda  Parameter  Estimates.  Table  6.5  lists  the  lambda 
parameter  estimations  of  the  “Normal”  distributed  data  samples  from  the  Monte  Carlo  sim¬ 
ulation  study.  The  first  row  contains  the  lambda  parameters  from  the  underlying  theoretical 
distribution.  The  next  nine  rows  contain  the  average  parameter  estimates  from  the  thirty 
samples  of  the  specified  set.  The  last  nine  rows  contain  the  variances  of  those  thirty  sample 


estimates. 


Fitment 

Method 


Theoretical 


Average  of 
Linear 
Moments 


Average  of 
Conventional 
Moments 


Average  of 

Alternate 

Moments 


Variance  of 
Linear 
Moments 


Variance  of 
Conventional 
Moments 


Variance  of 
Alternate 
Moments 


“Normal”  Distribution 


Sample 


oo 

0 

0.1975 

0.1349 

0.1349 

25 

-.2859 

0.1961 

0.1801 

0.2345 

50 

0.0269 

0.2368 

0.2386 

0.2032 

100 

-.0547 

0.2340 

0.1752 

0.1802 

25 

-.1466 

0.2427 

0.2055 

0.2544 

50 

0.0692 

0.2506 

0.2197 

0.2109 

100 

-.0462 

0.2515 

0.1884 

0.2047 

25 

0.0149 

0.2143 

0.2245 

0.1853 

50 

-.0303 

0.2288 

0.1974 

0.2061 

100 

-.0472 

0.2370 

0.1801 

0.1912 

0.7301 

0.0292 

0.0606 

0.0749 

0.3613 

0.0160 

0.0586 

0.0469 

0.1271 

0.0061 

0.0099 

0.0098 

0.8896 

0.0051 

0.0271 

0.0726 

0.3456 

0.0041 

0.0139 

0.0474 

0.1484 

0.0043 

0.0056 

0.0181 

0.7830 

0.0166 

0.0473 

0.0407 

0.5484 

0.0079 

0.0276 

0.0442 

0.1759 

0.0065 

0.0089 

0.0183 

Table  6.5  Lambda  Parameter  Estimates  of  “Normal”  Distributed  Samples 


Comparing  the  theoretical  values  to  the  estimated  sample  averages  shows  that  the 
method  of  linear  moments  usually  produces  the  most  accurate  estimates.  It  is  also  apparent 
that  an  increase  in  sample  size  does  not  necessarily  produce  an  increase  in  accuracy.  As 
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for  the  variances  of  the  sample  sets,  the  method  of  linear  moments  is  comparable  to  the 
method  of  conventional  moments  in  consistency  of  estimation.  Both  are  marginally  better 
than  alternate  moments.  In  all  cases,  an  increase  in  sample  size  produces  a  decrease  in 
variance,  suggesting  an  increase  in  the  consistency  of  parameter  estimation  for  all  methods, 
as  the  sample  size  increases. 

Table  6.6  lists  the  lambda  parameter  estimations  of  the  “Gamma”  distributed  data 
samples  from  the  Monte  Carlo  simulation  study.  In  this  case,  the  method  of  linear  moments 
again  produces  the  most  accurate  estimates.  The  averages  of  the  estimated  lambda  param¬ 
eters  are  closer  to  the  theoretical  values  than  are  the  averages  of  the  other  two  methods.  As 
for  consistency,  the  results  are  less  conclusive.  The  lowest  variances  are  scattered  randomly 
among  the  various  methods  and  sample  sizes.  Overall,  the  method  of  conventional  moments 
maintains  the  lowest  average  variance.  The  only  trend  is  that  the  variance  continues  to 
decrease  as  the  sample  size  increases. 
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Fitment 

Method 


“Gamma”  Distribution 


Sample 


■I 


Theoretical 

oo 

0 

.04134 

.00567 

.04046 

Average  of 

25 

-.1298 

0.0246 

-.0011 

0.0840 

Linear 

50 

-.1054 

0.0876 

0.0058 

0.1323 

Moments 

100 

-.0853 

0.0784 

0.0083 

0.0942 

Average  of 

25 

-.3565 

0.1927 

0.0132 

0.2916 

Conventional 

50 

-.2512 

0.1710 

0.0112 

0.2477 

Moments 

100 

-.2212 

0.1362 

0.0092 

0.1760 

Average  of 

25 

-.2003 

0.1361 

0.0158 

0.2065 

Alternate 

50 

-.1680 

0.1308 

0.0091 

0.1878 

Moments 

100 

-.0949 

0.0935 

0.0096 

0.1180 

Variance  of 

25 

0.0948 

0.0522 

0.0398 

Linear 

50 

0.0659 

0.0187 

0.0004 

0.0283 

Moments 

100 

0.0222 

0.0085 

0.0002 

0.0096 

Variance  of 

25 

0.0658 

0.0089 

0.0006 

0.0278 

Conventional 

50 

0.0570 

0.0057 

0.0002 

0.0239 

Moments 

100 

0.0343 

0.0045 

0.0001 

0.0099 

Variance  of 

25 

0.0782 

0.0263 

0.0010 

0.0428 

Alternate 

50 

0.0519 

0.0127 

0.0002 

0.0294 

Moments 

100 

0.0208 

0.0108 

0.0002 

0.0127 

Table  6.6  Lambda  Parameter  Estimates  of  “Gamma”  Distributed  Samples 

Table  6.7  lists  the  lambda  parameter  estimations  of  the  “Exponential"  distributed  data 
samples  from  the  Monte  Carlo  simulation  study.  In  this  case,  one  cannot  conclude  that  any 
method  is  substantially  more  accurate  t  han  the  others.  The  average  estimates  closest,  to  the 
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theoretical  values  are  spread  evenly  over  the  different  methods  and  sample  sizes.  However, 
these  results  indicate  the  alternate  moments  may  be  the  most  consistent  estimators.  The 
method  of  alternate  moments  has  the  lowest  variance  values  of  the  three  methods. 


Fitment 

Method 


Theoretical 


Average  of 
Linear 
Moments 


Average  of 
Conventional 
Moments 


Average  of 
Alternate 
Moments 


“Exponential” 


Sample 


Distribution 


HHI 


25  0.03S6 


50  0.0099 


-.0087 


25  0.1179 


50 


-.00163 


-.0794 


0.0277 


0.0170 


0.1522 


jXEfiaEXE 


100  0.1486 


-.0185 


25 


50 


100  -.0196 


0.0558 


0.0502 


iTtTiTiMmiMiill 


0.0426 


-.000009 


-.0064 


-.0013 


-.0002 


0.1756 


-.0001 


-.0001 


-.0012 


0.0000072 


0.000071 


-.00162 


-.0187 


0.0515 


0.0272 


0.1900 


0.1309 


0.0722 


0.0905 


0.0711 


0.0501 


Variance  of 

25 

0.0345 

0.0624 

0.0004 

0.0344 

Linear 

50 

0.0168 

0.0203 

0.00004 

0.0201 

Moments 

100 

0.0041 

0.0100 

0.000008 

0.0085 

Variance  of 

25 

0.2459 

0.0376 

0.4469 

0.0403 

Conventional 

50 

0.1000 

0.0113 

.00000003 

0.0257 

Moments 

100 

0.0654 

0.0080 

.00000002 

0.0128 

Variance  of 

25 

0.0246 

0.0354 

0.000084 

0.0365 

Alternate 

50 

0.0059 

0.0126 

0.0000037 

0.0151 

Moments 

100 

0.0024 

0.0060 

0.0000011 

0.0066 

Table  6.7  Lambda.  Parameter  Estimates  of  "Exponential”  Distributed  Samples 
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As  noted  for  the  previous  distributions,  the  variance  decreases  as  the  sample  size  in¬ 
creases,  for  all  three  methods.  The  lambda  parameter  values  for  the  theoretical  “Gamma” 
distribution  axe  about  three  times  smaller  than  the  parameter  values  of  the  theoretical  "Nor¬ 
mal”  distribution,  and  the  parameter  values  of  the  theoretical  “Exponential”  distribution  are 
about  twenty  times  smaller  than  the  parameters  of  the  “Gamma”  distribution.  Thus,  a  given 
change  in  a  lambda  parameter  is  more  likely  to  significantly  alter  the  appearance  of  an  esti¬ 
mated  “Exponential”  fit  than  an  estimated  “Normal”  fit,  even  though  the  change  may  not 
affect  the  K-S  values  for  goodness-of-fit  to  the  same  degree.  This  can  account  for  the  dis¬ 
crepancies  noted  between  the  visual  PDF’s  and  the  Kolmogorov-Smirnov  statistics  from  the 
comparable  CDF’s.  The  PDF’s  showed  that  the  number  of  poor  fits  increased  from  “Nor¬ 
mal”  to  “Gamma”  to  “Exponential”  shapes.  However,  the  averages  of  the  K-S  statistics 
decreased  along  the  same  line  of  comparison.  Even  though  the  K-S  statistics  were  small, 
the  difference  in  lambda  parameters  needed  to  change  a  “good  fit”  to  a  “poor  fit”  is  even 
smaller. 

6.3  Fitted  Distributions  vs.  Empirical  Distributions 

Until  now,  the  results  have  examined  the  three  methods’  abilities  to  closely  fit  a  distri¬ 
bution  to  the  underlying  theoretical  distribution.  The  methods  estimate  a  distribution  via 
sample  data,  rather  than  directly  from  the  theoretical  distribution.  This  intermediate  step  is 
the  source  of  the  “sampling  variability,”  or  “noise.”  We  can  now  also  examine  how  well  each 
of  the  methods  fit  a  distribution  to  the  empirical  CDF  of  the  individual  data  samples.  This 
ability  is  often  of  interest  to  a  modeler,  who  is  concerned  with  closely  fitting  a  distribution 
to  sample  data  when  the  underlying  distribution  of  the  population  is  unknown  and  is  used 
as  a  barometer  for  assessing  goodness-of-  fit. 

6.3. 1  Summarizing  Kolmogorov-  Smirnov  Statistics.  Table  6.8  summarizes  Kolmogorov  - 
Smirnov  (K-S)  statistics  collected  from  each  of  the  27  Monte  Carlo  experiments.  The  “Aver¬ 
age  K-S”  Column  lists  the  average  over  the  80  samples  of  the  average  difference  between  the 
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fitted  and  empirical  distribution  functions.  The  “Maximum  K-S”  Column  lists  the  average 
over  the  30  samples  of  the  maximum  difference  between  the  same  functions.  The  “Max  K-S 
Variance”  Column  lists  the  variance  of  the  30  values  which  form  the  average  in  the  previous 
column. 

Unlike  the  K-S  statistics  collected  from  the  theoretical  distributions,  these  statistics 
indicate  small  but  consistent  advantages  for  the  linear  moments  and  the  alternate  moments 
over  the  conventional  moments.  Regarding  the  variances  of  the  maximum  I\-S  values,  no 
method  shows  a  clear  advantage  for  “Normal”  distributed  samples,  but  both  the  linear 
and  alternate  methods  average  only  half  as  much  variance  as  the  conventional  moments  for 
“Gamma”  and  “Exponential”  distributed  sample  data.  With  both  consistently  lower  means 
and  variances  of  the  K-S  values,  the  methods  of  linear  and  alternate  moments  are  superior 
to  conventional  moments  for  estimating  and  fitting  distributions  to  sample  data. 
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Distribution 

Shape 

Sample 

Size 

Method 

of  Moments 

Average 

K-S 

Maximum 

K-S 

Max  K-S 

Variance 

linear 

0.02055 

0.09573 

.00052 

25 

conventional 

0.02146 

0.10001 

.00055 

alternate 

0.02042 

0.09788 

.00047 

linear 

0.01438 

0.06771 

.00031 

“Normal” 

50 

conventional 

0.01576 

0.07544 

.00061 

alternate 

0.01437 

0.07069 

.00036 

linear 

0.01001 

0.04640 

.00014 

100 

conventional 

0.01142 

0.05237 

.00025 

alternate 

0.01039 

0.05005 

.00024 

linear 

0.01720 

0.09493 

.00063 

25 

conventional 

0.01899 

0.11345 

.00096 

alternate 

0.01731 

0.10338 

.00074 

linear 

0.01198 

0.06632 

.00030 

“Gamma” 

50 

conventional 

0.01552 

0.08983 

.00079 

alternate 

0.01200 

0.06965 

.00034 

linear 

0.00855 

0.04647 

.00019 

100 

conventional 

0.01272 

0.07009 

.00040 

alternate 

0.00871 

0.04963 

.00024 

linear 

0.01841 

0.09955 

.00074 

25 

conventional 

0.01983 

0.13500 

.00138 

alternate 

0.01767 

0.11254 

.00130 

linear 

0.01196 

0.07062 

.00035 

“Exponential” 

50 

conventional 

0.01524 

0.10608 

.00110 

alternate 

0.01239 

0.07479 

.00035 

linear 

0.00855 

0.04811 

.00021 

100 

conventional 

0.01257 

0.08529 

.00083 

alternate 

0.00907 

0.05183 

.00040 
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Average  K-$  Values 
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Figure  6.16  Average  K-S  Values  of  fitted  CDF’s  vs.  empirical  CDF’s 

Figures  6.18  through  6.20  display  the  information  from  the  “Average  K-S,”  “Maximum 
K-S.”  and  “Max  K-S  Variance”  columns  of  Table  6.8,  respectively. 

Figure  6.18  clearly  shows  that  the  linear  and  alternate  moments  are  comparable  while 
the  conventional  moments  produce  slightly  higher  averages  of  the  K-S  values.  There  are  two 
downward  trends  depicted  in  the  chart.  The  average  K-S  values  decrease  as  the  sample  size 
increases,  and  as  the  skewness  and  kurtosis  of  the  theoretical  distribution  increases. 

Figure  6. 19  shows  even  more  disparity  between  the  conventional  moments  and  the  other 
two  methods.  Note  that  for  all  cases  of  “Normal’-  distributed  sample  data,  all  three  methods 
perform  equally  well.  For  “Gamma”  and  “Exponential"  distributed  sample  data,  the  K-S 
values  associated  with  teh  method  of  conventional  moments  become  incresingly  larger  than 
those  for  the  other  methods  as  the  skewness  and  kurtosis  increases.  Also,  for  the  first  time, 
linear  moments  show  a  small  but  consistent  improvement  over  alternate  moments.  Figure 
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Maximum  K-S  Values 


— ■ —  Linear  Moments 

— ° —  Conventional 
Moments 

— • —  Alternate  Moments 


Figure  6.17  Maximum  K-S  Values  of  fitted  CDF's  vs.  empirical  CDF’s 


6.19  also  shows  a  reversal  of  a  previous  trend:  the  K-S  values  increase  from  the  "Normal"’ 
to  the  "Gamma.”  and  from  the  “Gamma”  to  the  “Exponential”  sample  data,  for  a  constant 
sample  size. 

Figure  6.20  indicates  that  all  three  methods  perform  equally  well  for  all  cases  of  the 
“Normal”  distribution.  Alternate  moments  are  less  consistent  than  linear  moments  for  25- 
element  samples  but  are  equally  consisted  for  larger  samples.  Conventional  moments  show 
increasing  variance  for  “Gamma”  and  “Exponential"  distributions.  The  K-S  values  increase 
with  increasing  values  of  skewness  and  kurtosis. 

The  increasing  differences  in  the  K-S  values  also  suggest  that  the  method  of  conven¬ 
tional  moments  becomes  less  suitable  as  the  skewness  and  kurtosis  of  the  distribution  of 
sample  data  increase.  There  is  little  change  in  diffe.ences  between  K-S  values  as  the  sample 
size  changes,  as  depicted  by  the  nearly  parallel  lines  on  the  chart.  This  suggests  that,  while 
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Figure  6.  IS  Variances  of  Maximum  K-S  Values  of  fitted  CDF’s  vs.  empirical  CDF's 

the  methods  of  linear  and  alternate  moments  are  better  than  the  method  of  conventional 
moments  overall,  their  relative  superiority  is  less  for  small  data  samples.  Also,  this  does  not 
necessarily  indicate  an  advantage  in  fitting  distributions  to  the  underlying  distribution  of 
the  population. 

Table  6.S  and  Figures  6.18  through  6.20  provide  a  cursory  overall  look  at  the  objective 
measures  of  the  good ness-of-  fit  for  each  of  the  three  methods  of  fit.  The  mean  and  variance 
of  the  K-S  statistics  provide  a  feel  for  which  method  will  perform  well  in  the  long  run. 
and  will  provide  consistently -close  fits.  But,  information  on  which  method  will  most  oftei. 
provide  the  closest  fit  is  also  useful. 

6.3.2  Kohnogorov-Smirnov  Statistics  of  Individual  Data  Samples.  Tables  6.9  and 
6.10  list  the  number  of  times  that  each  method  of  fit  scored  the  lowest  K-S  value  among  each 
in  a  set  ol  thirty  samples.  I  he  tallies  are  divided  into  summaries  of  maximum  and  average 
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K-S  values.  Each  table  is  further  broken  down  by  the  underlying  distribution  and  size  of  the 
samples  in  the  set. 


Fitted  Distribution  vs.  Empirical  Distribution 

MAX  K-S  of 

MAX  K-S  of 

MAX  K-S  of 

Distribution  Shape 

Sample  Size 

Linear 

Conventional 

Alternate 

Method 

Method 

Method 

25 

14 

8 

S 

“Normal” 

50 

16 

5 

9 

100 

20 

4 

6 

25 

21 

4 

5 

“Gamma” 

50 

20 

2 

S 

100 

17 

3 

10 

25 

IS 

5 

7 

“Exponential” 

50 

20 

1 

9 

100 

14 

2 

14 

Table  6.9  Number  of  Lowest  MAX  K-S  Empirical  Differences  per  set  of  30  samples 


Table  6.9  shows  that  the  method  of  linear  moments  most  often  est  imates  the  best  fit  to 
the  distribution  of  the  sample  data  for  all  distribution  shapes.  The  linear  moments  show  clear 
superiority  over  both  of  the  other  methods.  However,  no  method  shows  consistent  increasing 
or  decreasing  performance  as  the  sample  size  changes  independent  of  the  sample  distribution. 
Most  notable  is  the  clear  advantage  that  the  method  of  linear  moments  possesses  in  light  of 
the  lack  of  similar  evidence  in  Table  6.3. 
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Fitted  Distribution  vs.  Empirical  Distribution 

AVG  K-S  of 

AVG  K-S  of 

AVG  K-S  of 

Distribution  Shape 

Sample  Size 

Linear 

Conventional 

Alternate 

Method 

Method 

Method 

25 

S 

12 

10 

“Normal” 

50 

8 

11 

11 

100 

19 

3 

8 

25 

12 

7 

11 

“Gamma” 

50 

14 

2 

14 

100 

14 

2 

14 

25 

9 

6 

15 

“Exponential” 

50 

14 

6 

10 

100 

16 

3 

11 

Table  6.10  Number  of  Lowest  AVG  K-S  Empirical  Differences  per  set  of  30  samples 


Table  6.10  shows  that  none  of  the  methods  are  better  estimators  of  the  distribution 
of  “Normal'’  distributed  sample  data  than  any  other.  However,  it  does  show  that  both  the 
method  of  linear  moments  and  the  method  of  alternate  moments  are  clearly  better  estimators 
of  empirical  distributions  than  the  method  of  conventional  moments.  Still,  there  is  no  change 
in  performance  which  corresponds  to  changes  in  the  sample  size  that  is  independent  of  the 
distribution  shape. 

Once  the  information  was  gleaned  from  the  data  collected  from  the  Monte  Carlo  ex¬ 
periment,  some  insight  was  gained  into  the  method  of  linear  moments.  These  findings  and 
the  conclusions  drawn  from  them  are  presented  in  Chapter  Seven. 
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VII.  Conclusions  and  Recommendations 


7.1  The  Generalized  Lambda  Distribution  and  Linear  Moments 

The  purpose  of  thesis  was  to  assess  the  usefulness  of  the  method  of  L-moments  applied 
to  the  GLD,  especially  in  relation  to  other,  more  well-established  methods  for  parameter 
estimation.  The  specific  goal  was  to  determine  whether  any  one  method  is  superior  in 
fitting  the  generalized  lambda  distribution  to  sample  data  when  the  underlying  distribution 
is  unknown  and  the  sample  size  is  relatively  small. 

The  results  in  Chapter  Six  suggest  that  the  methods  of  linear  and  alternate  moments 
perform  about  equally  well  and  are  both  superior  to  the  conventional  method  of  moments 
for  fitting  the  estimated  GLD  function  to  an  underlying  theoretical  distribution.  They 
also  indicate  that  the  methods  of  linear  and  alternate  moments  perform  equally  well.  This 
superiority  is  evident  in  the  low  averages  and  variances  of  the  Kolmogorov-Smirnov  statistics, 
the  accurate  averages  and  low  variances  of  the  estimated  lambda  parameters,  and  the  high 
percentages  of  properly-shaped  PDF’s  generated  from  the  sets  of  sample  data.  Thus,  linear 
moments  can  be  considered  an  improvement  on  the  methods  of  conventional  and  alternate 
moments. 

The  results  in  Chapter  Six  also  suggest  that  the  method  of  linear  moments  is  better 
than  the  other  two  methods  tested  for  fitting  the  GLD  to  the  empirical  distributions  of 
sample  data.  This  superiority  is  evident  in  the  low  averages  and  variances  of  the  Kolmogorov- 
Smirnov  statistics  collected  from  comparisons  of  the  fitted  distributions  and  the  empirical 
distributions  of  sample  data. 

In  addition,  the  method  of  linear  moments  showed  the  greatest  flexibility  in  fitting 
the  GLD  to  sample  data  from  different  distributions.  The  method  of  L-moments  estimated 
the  GLD  from  "Gamma"  distributed  sample  data  as  well  as  it  did  from  normal-distributed 
sample  data,  and  there  was  only  moderate  degradation  in  the  estimates  of  "Exponential” 
distributee!  samples.  The  method  of  alternate  moments  showed  increasing  difficulty  with 


both  the  '‘Gamma”  and  “Exponential”  samples.  The  conventional  method  of  moments  had 
difficulty  when  dealing  with  “Gamma”  distributed  and  “Exponential”  distributed  sample 
data.  It  appears  that  the  asymmetry  and  convexity,  as  measured  by  the  skewness  and 
kurtosis,  respectively,  of  sample  data  can  also  have  a  negative  effect  on  the  good  ness-of- fit 
for  some  methods  of  estimation,  such  as  the  conventional  method  of  moments.  The  results 
of  this  experiment  show  that  linear  moments  and  alternate  moments  are  much  less  affected 
by  the  asymmetry  and  convexity  of  the  empirical  distributions  of  sample  data. 

In  general,  the  method  of  linear  moments  shows  an  increasing  advantage  as  the  asym¬ 
metry  and  convexity  of  the  distribution  of  sample  data  increase,  and  a  constant  but  small 
advantage  over  the  other  methods  tested  as  the  sample  size  decreases.  However,  these  advan¬ 
tages  are  gained  primarily  when  fitting  estimated  distributions  to  the  empirical  distribution 
of  an  individual  sample.  The  advantage  was  less  evident  when  using  linear  moments  to  es¬ 
timate  the  underlying  theoretical  distribution  from  sample  data  than  when  using  the  other 
two  methods. 

I  recommend  that  those  involved  in  modeling  and  simulation  use  linear  moments  as 
the  preferred  method  when  fitting  the  Generalized  Lambda  Distribution  to  sample  data.  In 
addition  to  my  findings  and  conclusions,  I  have  recommendations  for  further  research  in  the 
following  areas. 

1.2  GLD  Software  Package,  Version  II,  in  C++ 

Hsu  (1991)  developed  a  complete  software  package  using  the  GLD  and  various  esti¬ 
mating  methods  to  fit  distributions  to  sample  data  and  to  derive  lambda  parameters  from 
various  known  statistics.  The  software  package  is  written  in  C  -f  +.  a  state-of-the-art. 
object  -oriented  computer  language.  The  program  is  menu-driven,  user-interactive,  and 
employs  real  time  color  graphics.  The  underlying  algorithms  of  the  program  employ  the 
conventional  method  of  moments  and  the  alternate  method  of  moments,  as  well  as  oth¬ 
ers.  L-moments  should  be  incorporated  into  a  new  version  of  the  software  to  enhance  its 
capabilities. 


7 .3  Using  Powell’s  Algorithm  with  Hosking’s  Lambda  Distribution 

As  discussed  in  Chapter  Four,  Powell’s  Algorithm  is  the  numerical  search  routine  of 
choice  for  this  study  and  for  the  GLD.  However,  it  is  still  not  completely  proven  to  be  reliable 
when  searching  back  and  forth  between  the  positive  and  negative  regions.  Hosking  (August 
1993)  encountered  the  same  problem  using  Newton’s  method  and,  instead  of  revising  the 
search  routine  as  Mykytka  has  done,  he  revised  the  lambda  distribution  quantile  function. 
Hosking’s  Lambda  distribution  could  be  used  with  Powell’s  Algorithm  to  create  a  more 
reliable,  faster  search  routine.  This  combination  would  permit  the  search  routine  to  begin 
at  the  origin,  proceed  with  the  search  in  the  right  direction  the  first  time,  and  revisit  the 
origin  if  needed,  avoiding  the  forbidden  zones  during  transitions  between  regions. 

7.4  Expanding  the  GLD’s  Range  of  Approximate  Distribution  Shapes 

As  Hosking  noted,  the  sample  L-moments  exist  if  the  sample  has  a  finite  mean.  This 
does  not  always  hold  true  for  conventional  moments.  Thus,  it  is  possible  that  the  method  of 
linear  moments  will  permit  modelers  to  compute  statistics  from  samples  that  previously  had 
no  existing  conventional  moments.  Given  the  sample  L-moments.  GLD  parameters  could  be 
calculated,  and  a  distribution  could  be  fit.  These  fitted  distributions  may  extend  into  the 
regions  of  skewness  and  kurtosis  previously  unusable  by  the  GLD,  as  depicted  in  Figure  2.2. 


Appendix  A.  Comparisons  of  L-Moment  Equations 

In  this  Appendix,  I  present  an  in-depth  comparison  of  the  similarities  and  differences 
of  the  equations  for  the  four  L-moments,  Ai,A2,  A3,  and  A4  from  the  four  parameters  the 
Generalized  Lambda  Distribution  and  the  five  parameters  of  Hosking’s  Lambda  Distribution. 

I  begin  by  comparing  the  quantile  function  of  Hosking’s  Lambda  Distribution  to  the 
percentile  function  of  the  Generalized  Lambda  Distribution  to  determine  the  conversion 
factors  between  the  two  functions. 

Starting  with  Hosking’s  quantile  function  (Hosking,  1986:  83): 

x(F)  =  t  +  aF‘}  F)s,  0  <  F  <1.  (A.l) 

By  setting  a  =  7,  Hosking  converts  this  five  parameter  function  to  a  four  parameter 
function  (Hosking,  1993:  1): 

x{F)  =  t  +  aFp-  a(l  -  F)s,  0  <  F  <  1.  (A.2) 

I  now  introduce  a  new  term,  e,  which  equals  the  inverse  of  Hosking’s  a  term: 

a  =  ~>  (A-3) 

Fp  (1  -  F)s 

*(F)  =  ^  +  —  — ■^-L,  0  <  F  <  1.  (A. 4) 

Combining  the  two  terms  into  a  single  term  with  a  common  denominator: 


A-l 


x(F)  =  (  + 


F*  -  (1  -  F)s 


(A.5) 


,  0  <  F  <  1. 

£ 

Replacing  Hosking’s  variable  names  with  terminology  used  by  Ramberg  for  the  GLD, 
F  =  p, x(F)  —  R(p),(  =  Ai,e  =  Aj,/3  =  A3, 5  =  A4,  results  in  the  percentile  function  of  the 
GLD  (Bergevin,  1993:  7): 


«(p)  =  A ,  +  ***’  (1  P>>',  0<p<l.  (A.6) 

Aj 

Therefore,  the  two  operations  necessary  to  convert  Hosking’s  Lambda  Distribution  to 
Ramberg’s  Generalized  Lambda  Distribution  are  to  set  the  second  and  fourth  parameters, 
a  and  7  equal  to  each  other;  and  then  to  replace  them  with  their  inverse,  e.  These  two 
operations  will  be  used  to  show  the  similarity  between  L-moment  equations  of  Hosking’s 
Lambda  Distribution  and  the  GLD. 

Beginning  with  Hosking’s  equation  for  the  first  L-moment  (Hosking,  1986:  83): 


Ai  =  *  + 


a 


1 +/3  1  +  6 

Cross-multiplying  for  a  common  denominator: 


Ai=*  + 


a(l  +  6) 


7(1  +  P) 


(1+/?)(1  +f)  (l+f)(l  +  0)’ 

Combining  the  terms  with  a  common  denominator: 


Ai=e  + 


<*(!  +  *)  —  7(!+fl) 
(l  +  ^)(l+/3)  ' 


Replacing  a  and  7  with  e: 


(A.7) 


(A.8) 


(A.9) 


A-2 


(A10) 


(1  +  6)  —  (1  -f 

~U+S)(l+0) 


Simplifying  the  numerator: 


Ai  =  £  + 


6-0 


(All) 


e(l  +  6)(l  +  0Y 

Replacing  Hosking’s  terms  with  Ramberg’s  terms  yields  the  equation  for  Ax  derived 
by  Bergevin  (Bergevin,  1993:  64): 


Ax  —  Ax  + 


A<  —  A3 


(A. 12) 


Aj(A3  +  1)(A<  +  1) 

Performing  the  same  operations  on  Hosking’s  equation  for  the  second  L-moment  (Hosk- 
ing,  1986:  83): 


a0  7  6 

(l+0)(2  +  0)  +  (l  +  6)(2  +  8Y 


Cross-multiplying  for  a  common  denominator: 


(A-13) 


afl(l+f)(2  +  f)  7S(  l  +  fl)(2+0) 

(1  +  0)(2  +  /})(!  +  «)(2  +  *)  (!  +  *)(2  +  «)(!  +  m  +  «■ 


(A.14) 


Combining  the  terms  with  a  common  denominator: 


a0(l  +  S)(2  +  5)  +  7tf(l  +  0)(2  +  0) 
1  (l  +  «)(2  +  «)(l  +  0)(2+/?) 


(A.  15) 


Replacing  a  and  7  with  e: 


0(l  +  6)(2  +  6)  +  6(l+0)(2  +  0) 
1  e(l  +  5)(2  +  5)(1  +  0)(2  +  0) 


A-3 


Replacing  Hosking’s  terms  with  Ramberg’s  terms  yields  the  equation  for  A]  derived 
by  Bergevin  (Bergevin,  1993:  65): 

A  —  +  A4(A3  +  1)(^3  +  2) 

A2(A3  -f  1)(A3  -f  2)(A4  +  1)(A4  -+■  2) 

Performing  these  operations  on  Hosking’s  equation  for  the  third  L-moment  (Hosking, 
1986:  83): 


~  1) _ ~  1) 

(l-f/?)(2  +  /3)(3  +  /9)  (l  +  $)(2  +  5)(3  +  5)' 


Cross-multiplying  for  a  common  denominator: 


(A.18) 


a0(0-l)(l+S)(3+S)(3+6) _ yS(S~l)(l+/3)(3+0)(3+0) 

(i+/3)(a+fl)(s+/j)(i+<)(a+tf)(3+i)  (i+«)(a+«)(3+*)(i+/J)(a+0)(3+/8)’ 


(A.19) 


Combining  the  terms  with  a  common  denominator: 


1X3  ~  (l+«)(a+4)(3+«)(l+/3)(a+/3)(3+/9) 


(A. 20) 


Replacing  a  and  7  with  e: 


_  0(0-l)(l+6)(3+6)(3+f)-6(S-l)(l+0)(2+0)(3+0) 

«(i+tf)(a+i)(3+i)(i+/3)(a+ja)(3+/3) 


Simplifying  the  numerator: 


(A-21) 


A-4 


r 


A3  = 


(A. 22) 


Replacing  Hosking’s  terms  with  Ramberg’s  terms  yields  the  equation  for  Aa  derived 
by  Bergevin  (Bergevin,  1993:  66): 


A  _  (Al-Aa)(A44-X)(A4+a)(A4-H)-(A»-A4)(Aa-H  )(A»+2)(  Aa+3) 

*  Aj(Aj+l)(Aj+a)(Ai+3)(A4+l)(A4+J)(A4+3) 


(A. 23) 


Finally,  performing  these  operations  on  Hosking’s  equation  for  the  fourth  L-moment 
(Hosking,  1986:  83): 


~  (l+/3)(2+^)(3+j9)(4+/J)  ^  (l+*)(2+«)(3+S)(4+«)‘ 


Cross-multiplying  for  a  common  denominator: 


(A.24) 


1  isy-iKt-m+m+m+m+P) 

Ml+«)(2+tf)(3+«)(4+«)(l+0)(2+/9)(3+0)(4+0)- 


(A. 25) 


Combining  the  terms  with  a  common  denominator: 


A  _  c,0(0-l)(0-2)(l+6)(2+6)(3+6)(4+6)+‘l6l6-l)(6-2)(l+0)(2+0)(3+0)(i+0) 

V4  (l+«)(2+i)(3+i)(4+«)(l-H9)(2+0)(3+0)(4+0) 


(A. 26) 


Replacing  a  and  7  with  e: 


A-5 


Simplifying  the  numerator: 


*  _  ti9a-3fla+afl)(l+<)(a+M+m+<)+(g»-«3-a<)(l+fl)(a+/9)(3+fl)(4+fl) 

A4  “  e(l+«)(a+«)(3+«)(4+i)(l+^)(a+/3)(3+/3)(4+/9) 


(A. 28) 


Replacing  Hosking’s  terms  with  Ramberg’s  terms  yields  the  equation  for  A<  derived 
by  Bergevin  (Bergevin,  1993:  69): 


(A3-3A3+aAa)(A4+l)(A«+a)(A4+3)(A«+4)+(A»-3A3+aA«)(Ai+l)(Aa+a)(Ai+3)(AJ+4) 

Aa(As+l)(As+a)(A3+3)(A,+4)(A«+x)(A4+a)(A4+3)(A4+4) 


(A.29) 


Appendix  B.  Experiment  Results:  Normal  Distribution ,  Linear  Moments 

Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

-0.32920000 

0.33291600 

-0.08192267 

0.03075791 

Moments  of 

2nd 

0.46495600 

0.69995600 

0.56678553 

0.00435006 

Sample  Data 

3rd 

-0.14980300 

0.16846000 

0.02440728 

0.00632024 

4th 

-0.01568870 

0.20792600 

0.09713558 

0.00315052 

-2.07855000 

1.39302000 

-0.28595313 

0.73005042 

Lambda 

-0.22071600 

0.49963800 

0.19608461 

0.02915614 

Parameters 

-0.08825540 

1.00248000 

0.18008392 

0.06060880 

B 

-0.12274100 

1.11750000 

0.23446677 

0.07488433 

Theoretical 

MIN 

0.00000211 

0.00132966 

0.00015583 

0.00000008 

K-S 

AVG 

0.01000620 

0.05861190 

0.03572510 

0.00016930 

Statistics 

MAX 

0.02489270 

0.18884700 

0.10145870 

0.00149412 

Empirical 

MIN 

0.00000000 

0.00024583 

0.00002845 

0.00000000 

K-S 

AVG 

0.01290350 

0.02903130 

0.02055428 

0.00002040 

Statistics 

MAX 

0.06148730 

0.14373100 

0.09572504 

0.00051787 

Experiment  Number:  1A 

Samples:  30  samples  of  25  elements  each 

Original  function:  lambda  approximation  of  a  normal  function 
Method  of  fitment:  method  of  linear  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

-0.24640400 

0.34302800 

-0.04351340 

0.01692519 

Moments  of 

2nd 

0.46380700 

0.65934200 

0.56482997 

0.00248374 

Sample  Data 

3rd 

-0.12463400 

0.08293850 

-0.00737903 

0.00265956 

4th 

-0.00161649 

0.18735300 

0.10071303 

0.00163888 

-1.49791000 

1.53929000 

0.02692249 

0.36133215 

Lambda 

-0.12565200 

0.60877900 

0.23682851 

0.01601138 

Parameters 

-0.05730220 

1.12361000 

0.05862897 

H 

-0.05661900 

0.96255000 

0.20318173 

0.04691401 

Theoretical 

MIN 

0.00000482 

0.00128289 

0.00016679 

0.00000010 

K-S 

AVG 

0.00678552 

0.05570450 

0.02451674 

0.00010957 

Statistics 

MAX 

0.02189700 

0.13389000 

0.06728391 

0.00085707 

Empirical 

MIN 

0.00000000 

0.00007557 

0.00001988 

0.00000000 

K-S 

AVG 

0.01005300 

0.02212740 

0.01435771 

0.00001019 

Statistics 

MAX 

0.04356000 

0.12729100 

0.06770739 

0.00031097 

Experiment  Number:  IB 

Samples:  30  samples  of  50  elements  each 

Original  function:  lambda  approximation  of  a  normal  function 

Method  of  fitment:  method  of  linear  moments 
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Figure  B.2  Plots  of  50-element  samples  from  an  approximate  normal  distribution  fitted 
with  linear  moments 


Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

-0.21766800 

0.23810900 

-0.05008467 

0.01089713 

Moments  of 

2nd 

0.50551300 

0.61714400 

0.55958817 

0.00074421 

Sample  Data 

3rd 

-0.09907760 

0.08210260 

-0.00073607 

0.00193096 

4th 

0.06930190 

0.16791000 

0.11116174 

0.00041959 

D 

-0.65020500 

0.74872200 

-0.05474650 

0.12705271 

Lambda 

-0.00664537 

0.38615300 

0.23400949 

0.00608533 

Parameters 

-0.00288072 

0.49463600 

0.17515935 

0.00994712 

H 

-0.00406039 

0.40454100 

0.18023218 

0.00977001 

Theoretical 

MIN 

0.00000004 

0.00123381 

0.00021441 

0.00000016 

K-S 

AVG 

0.00733220 

0.03906030 

0.01885646 

0.00007125 

Statistics 

MAX 

0.01486220 

0.08808770 

0.05100268 

0.00051479 

Empirical 

MIN 

0.00000000 

0.00011953 

0.00002576 

0.00000000 

K-S 

AVG 

0.00657308 

0.01709660 

0.01000614 

0.00000579 

Statistics 

MAX 

0.02922960 

0.08260010 

0.04640222 

0.00014395 

Experiment  Number:  1C 

Samples:  30  samples  of  100  elements  each 

Original  function:  lambda  approximation  of  a  normal  function 

Method  of  fitment:  method  of  linear  moments 
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Figure  B.3  Plots  of  100-element  samples  from  an  approximate  normal  distribution  fitted 
with  linear  moments 
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Appendix  C.  Experiment  Results:  Gamma  Distribution,  Linear  Moments 

Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

II 

Statistic  j 

Minimum 

Maximum 

Average 

Variance 

1st 

0.48061600 

1.171880C0 

0.74208173 

0.03016900 

Moments  of 

2nd 

0.38186000 

0.67193700 

0.51891087 

0.00718194 

Sample  Data 

3rd 

0.07615490 

0.37085800 

0.23640660 

0.00690586 

4th 

0.02031400 

0.29404600 

0.14745152 

0.00407120 

-0.62057600 

0.49394800 

-0.12984145 

0.09483345 

Lambda 

-0.65561600 

0.31523900 

0.02458554 

0.05222803 

Parameters 

-0.10849900 

0.05369720 

-0.00112102 

0.00122820 

H 

-0.32500000 

0.51999700 

0.08401853 

0.03978563 

Theoretical 

MIN 

0.00000067 

0.00088185 

0.00011831 

0.00000004 

K-S 

AVG 

0.01115250 

0.04637010 

0.02790449 

0.00010405 

Statistics 

MAX 

0.02728510 

0.16391300 

0.09577868 

0.00132345 

Empirical 

MIN 

0.00000000 

0.00022676 

0.00002922 

0.00000000 

K-S 

AVG 

0.01051120 

0.02525230 

0.01720370 

0.00001793 

Statistics 

MAX 

0.05963730 

0.15913600 

0.09493287 

0.00063476 

Experiment  Number:  2A 

Samples:  30  samples  of  25  elements  each 

Original  function:  lambda  approximation  of  a  gamma  function 
Method  of  fitment:  method  of  linear  moments 
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Figure  C.l  Plots  of  25-element  samples  from  an  approximate  gamma  distribution  fitted 
with  linear  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.55609200 

1.23099000 

0.76728553 

0.01912877 

Moments  of 

2nd 

0.40011800 

0.73080100 

0.51610747 

0.00496379 

Sample  Data 

3rd 

0.10441400 

0.31780000 

0.20907797 

0.00304550 

4th 

0.04483840 

0.22731200 

0.13154492 

0.00207545 

D 

-0.59335000 

0.36391500 

-0.10537867 

0.06585115 

Lambda 

1 

-0.35408300 

0.30237500 

0.08758405 

0.01873872 

Parameters 

-0.06458690 

0.04024890 

0.00584606 

0.00037417 

H 

-0.17902700 

0.49252300 

0.13232465 

0.02834047 

Theoretical 

MIN 

0.00000033 

0.00088185 

0.00015311 

0.00000006 

K-S 

AVG 

0.00420067 

0.05758760 

0.02040218 

0.00009339 

Statistics 

MAX 

0.02140450 

0.12633600 

0.06529886 

0.00076624 

Empirical 

MIN 

0.00000000 

0.00011152 

0.00001303 

0.00000000 

K-S 

AVG 

0.00687171 

0.01760340 

0.01197734 

0.00000882 

Statistics 

MAX 

0.04019940 

0.11113000 

0.06632050 

0.00030197 

Experiment  Number:  2B 

Samples:  30  samples  of  50  elements  each 

Original  function:  lambda  approximation  of  a  gamma  function 

Method  of  fitment:  method  of  linear  moments 
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Figure  C.2  Plots  of  50-element  samples  from  an  approximate  gamma  distribution  fitted 
with  linear  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.57285600 

1.09224000 

0.76092617 

0.01117079 

Moments  of 

2nd 

0.41755600 

0.66484400 

0.51482683 

0.00213806 

Sample  Data 

3rd 

0.12744100 

0.29910900 

0.22375943 

0.00187228 

4th 

0.07923440 

0.22313900 

0.14053823 

0.00079948 

D 

-0.47898200 

0.17442600 

-0.08526832 

0.02220021 

Lambda 

-0.25092900 

0.21989700 

0.07842475 

0.00847902 

Parameters 

-0.03406480 

0.03163950 

0.00830060 

0.00017585 

H 

-0.16210500 

0.28390000 

0.09417172 

0.00964847 

Theoretical 

MIN 

0.00000111 

0.00088185 

0.00014947 

0.00000006 

K-S 

AVG 

0.00632982 

0.03989130 

0.01549239 

0.00005874 

Statistics 

MAX 

0.01531930 

0.08619750 

0.05020063 

0.00043966 

Empirical 

MIN 

0.00000000 

0.00016237 

0.00002068 

0.00000000 

K-S 

AVG 

0.00517  n3 

0.01353660 

0.00855233 

0.00000408 

Statistics 

MAX 

0.02793380 

0.08175310 

0.04647265 

0.00018876 

Experiment  Number:  2C 

Samples:  30  samples  of  100  elements  each 

Original  function:  lambda  approximation  of  a  gamma  function 

Method  of  fitment:  method  of  linear  moments 
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Figure  C.3  Plots  of  100-element  samples  from  an  approximate  gamma  distribution  fitted 
with  linear  moments 
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Appendix  D.  Experiment  Results:  Exponential  Distribution ,  Linear 

Moments 

Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.67087600 

1.36818000 

0.93589840 

0.02967292 

Moments  of 

2nd 

0.32504700 

0.68764400 

0.48435410 

0.00919645 

Sample  Data 

3rd 

0.19025100 

0.46263500 

0.33083303 

0.00652701 

4th 

0.02826610 

0.32878800 

0.16516392 

0.00527456 

-0.24390200 

0.54126400 

0.03855568 

0.03454463 

Lambda 

-0.78781700 

0.22410700 

-0.07935519 

0.06237191 

Parameters 

1 

-0.08832810 

0.01016820 

-0.00644237 

0.00040783 

H 

-0.37972500 

0.36723100 

-0.01870103 

0.03442562 

Theoretical 

0.00000000 

0.00127325 

0.00022416 

0.00000012 

K-S 

AVG 

0.00905325 

0.04554390 

0.02479719 

0.00009789 

Statistics 

MAX 

0.03506210 

0.16055200 

0.09236931 

0.00127595 

Empirical 

MIN 

0.00000000 

0.00021511 

0.00006333 

0.00000000 

K-S 

AVG 

0.00968736 

0.04757100 

0.01840903 

0.00006827 

Statistics 

MAX 

0.05938320 

0.18146000 

0.09954951 

0.00074020 

Experiment  Number:  3 A 

Samples:  30  samples  of  25  elements  each 

Original  function:  lambda  approximation  of  an  exponential  function 
Method  of  fitment:  method  of  linear  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.75096200 

1.44070000 

0.95519720 

0.01969888 

Moments  of 

2nd 

0.35598200 

0.75015500 

0.48163140 

0.00651215 

Sample  Data 

3rd 

0.22137500 

0.42617500 

0.30446793 

0.00300091 

4th 

0.04726550 

0.24158300 

0.14203869 

0.00257635 

D 

-0.20071400 

0.33470900 

0.00991586 

0.01683758 

Lambda 

1 

-0.39034500 

0.26419000 

0.02766978 

0.02028058 

Parameters 

-0.03202840 

0.00766124 

-0.00126660 

0.00004343 

H 

-0.20218900 

0.32630700 

0.05152388 

0.02006486 

Theoretical 

MIN 

0.00000074 

0.00187586 

0.00034786 

0.00000026 

K-S 

AVG 

0.00366195 

0.06010810 

0.01895440 

0.00010956 

Statistics 

MAX 

0.02063870 

0.11822400 

0.06403366 

0.00063196 

Empirical 

MIN 

0.00000000 

0.00028062 

0.00002655 

0.00000000 

K-S 

AVG 

0.00600925 

0.01978000 

0.01196005 

0.00001157 

Statistics 

MAX 

0.03931390 

0.10229500 

0.07062143 

0.00035156 

Experiment  Number:  3B 

Samples:  30  samples  of  50  elements  each 

Original  function:  lambda  approximation  of  an  exponential  function 
Method  of  fitment:  method  of  linear  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.76118800 

1.29035000 

0.94777557 

0.01100065 

Moments  of 

2nd 

0.36137600 

0.66616900 

0.48157510 

0.00302355 

Sample  Data 

3rd 

0.24242200 

0.40198900 

0.32050903 

0.00173883 

4th 

0.09137160 

0.24539800 

0.15242416 

0.00109253 

D 

-0.11582100 

0.15770800 

-0.00874052 

0.00410422 

Lambda 

. 

-0.30981200 

0.17745500 

0.01700511 

0.00996706 

Parameters 

-0.00982581 

0.00567225 

-0.00018382 

0.00000784 

H 

-0.19996100 

0.22428400 

0.02720114 

0.00850612 

Theoretical 

MIN 

0.00000037 

0.00635650 

0.00041018 

0.00000134 

K-S 

AVG 

0.00609105 

0.04133410 

0.01449844 

0.00006345 

Statistics 

MAX 

0.01750830 

0.08462420 

0.04547370 

0.00038172 

Empirical 

MIN 

0.00000000 

0.00012106 

0.00003004 

0.00000000 

K-S 

AVG 

0.00473884 

0.01318650 

0.00855294 

0.00000501 

Statistics 

MAX 

0.02954890 

0.08860690 

0.04811859 

0.00021221 

Experiment  Number:  3C 

Samples:  30  samples  of  100  elements  each 

Original  function:  lambda  approximation  of  an  exponential  function 
Method  of  fitment:  method  of  linear  moments 
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Appendix  E.  Experiment  Results:  Normal  Distribution,  Conventional 

Moments 

Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

-0.32920000 

0.33291600 

-0.08192267 

0.03075791 

Moments  of 

2nd 

0.68745300 

1.54415000 

1.03898053 

0.05250603 

Sample  Data 

3rd 

-0.98539500 

0.63038300 

-0.05117874 

0.16026090 

4th 

1.83407000 

4.16280000 

2.67890867 

0.30904123 

D 

-1.83070000 

1.52185000 

-0.14659537 

0.88960911 

Lambda 

1 

0.07862730 

0.38893800 

0.24267934 

0.00508416 

Parameters 

0.00000526 

0.54493000 

0.20546724 

0.02714083 

H 

0.00562861 

1.16183000 

0.25439328 

0.07261522 

Theoretical 

MIN 

0.00000242 

0.00132966 

0.00022295 

0.00000013 

K-S 

AVG 

0.01538140 

0.05529040 

0.03547728 

0.00015160 

Statistics 

MAX 

0.04246670 

0.16926800 

0.09836652 

0.00104047 

Empirical 

MIN 

0.00000000 

0.00008005 

0.00001482 

0.00000000 

K-S 

AVG 

0.01102580 

0.03059050 

0.02146268 

0.00003285 

Statistics 

MAX 

0.05243800 

0.14790400 

0.10000633 

0.00054649 

Experiment  Number:  4 A 

Samples:  30  samples  of  25  elements  each 

Original  function:  lambda  approximation  of  a  normal  function 
Method  of  fitment:  conventional  method  of  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

-0.24640400 

0.34302800 

-0.03126540 

0.02190964 

Moments  of 

2nd 

0.67163200 

1.32877000 

1.01719943 

0.03144692 

Sample  Data 

3rd 

-0.70625700 

0.28223700 

-0.12639136 

0.06207302 

4th 

2.07060000 

3.46496000 

2.70478833 

0.12615150 

-1.20415000 

1.24473000 

0.06923864 

0.34558201 

Lambda 

0.10878000 

0.38976400 

0.25059593 

0.00405041 

Parameters 

1 

0.06610440 

0.45249000 

0.21968641 

0.01388314 

H 

0.03763780 

1.03470000 

0.21087366 

0.04743812 

Theoretical 

MIN 

0.00000043 

0.00132966 

0.00017865 

0.00000012 

K-S 

AVG 

0.01055100 

0.05653010 

0.02521661 

0.00014575 

Statistics 

MAX 

0.02042070 

0.12527700 

0.06705405 

0.00100233 

Empirical 

MIN 

0.00000000 

0.00022214 

0.00003048 

0.00000000 

K-S 

AVG 

0.01003980 

0.03600580 

0.01576388 

0.00002948 

Statistics 

MAX 

0.04035370 

0.15387900 

0.07544494 

0.00061213 

Experiment  Number:  4B 

Samples:  30  samples  of  50  elements  each 

Original  function:  lambda  approximation  of  a  normal  function 

Method  of  fitment:  conventional  method  of  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

-0.21766800 

0.23810900 

-0.04653627 

0.01058862 

Moments  of 

2nd 

0.79556500 

1.19861000 

0.99209797 

0.00903663 

Sample  Data 

3rd 

-0.57158700 

0.23436100 

-0.04420776 

0.03683652 

4th 

2.06409000 

3.64464000 

2.76636533 

0.09748634 

D 

-1.17677000 

0.91700300 

-0.04620887 

0.14842199 

Lambda 

0.09009060 

0.39928400 

0.25147062 

0.00426265 

Parameters 

■ 

0.04436510 

0.38893000 

0.18841034 

0.00562870 

H 

0.05461550 

0.70557900 

0.20469460 

0.01811967 

Theoretical 

MIN 

0.00000005 

0.00132745 

0.00019018 

0.00000015 

K-S 

AVG 

0.00780327 

0.03940710 

0.01836428 

0.00007217 

Statistics 

MAX 

0.01603560 

0.08966780 

0.04879070 

0.00048041 

Empirical 

MIN 

0.00000000 

0.00008975 

0.00002126 

0.00000000 

K-S 

AVG 

0.00642899 

0.02660000 

0.01141977 

0.00001761 

Statistics 

MAX 

0.02808660 

0.08697620 

0.05236926 

0.00024790 

Experiment  Number:  4C 

Samples:  30  samples  of  100  elements  each 

Original  function:  lambda  approximation  of  a  normal  function 

Method  of  fitment:  conventional  method  of  moments 
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Figure  E.3  Plots  of  100-element  samples  from  an  approximate  normal  distribution  fitted 
with  conventional  moments 
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Appendix  F.  Experiment  Results:  Gamma  Distribution ,  Conventional 

Moments 

Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.48061600 

1.17188000 

0.74208173 

0.03016900 

Moments  of 

2nd 

0.44652400 

1.83663000 

0.93963863 

0.11671945 

Sample  Data 

3rd 

0.27005400 

2.16257000 

0.99939517 

0.27040342 

4th 

1.95383000 

8.77330000 

3.77792533 

2.94413918 

D 

-0.77041100 

0.15703100 

-0.35654108 

0.06580892 

Lambda 

0.00431302 

0.38850300 

0.19268309 

0.00889336 

Parameters 

D 

0.00000020 

0.09125090 

0.01315165 

0.00059144 

H 

0.00393935 

0.66545900 

0.29161819 

0.02784793 

Theoretical 

MIN 

0.00000073 

0.00088285 

0.00023834 

0.00000007 

K-S 

AVG 

0.01365960 

0.05092160 

0.02837313 

0.00011045 

Statistics 

MAX 

0.03297290 

0.16598500 

0.09354674 

0.00103906 

Empirical 

MIN 

0.00000000 

0.00039474 

0.00002796 

0.00000001 

K-S 

AVG 

0.00998735 

0.04265360 

0.01898630 

0.00005508 

Statistics 

MAX 

0.05278450 

0.18413800 

0.11344866 

0.00095563 

Experiment  Number:  5A 

Samples:  30  samples  of  25  elements  each 

Original  function:  lambda  approximation  of  a  gamma  function 
Method  of  fitment:  conventional  method  of  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.55609200 

1.23099000 

0.76728553 

0.01912877 

Moments  of 

2nd 

0.51725800 

1.95520000 

0.08046629 

Sample  Data 

3rd 

0.43535500 

1.50274000 

0.97441073 

0.10478642 

4th 

2.22859000 

5.88553000 

3.78924267 

1.14210294 

D 

-0.63880000 

0.18232200 

-0.25121066 

0.05703032 

Lambda 

0.04615670 

0.34312800 

0.17102693 

0.00565214 

Parameters 

1 

0.00000032 

0.06328690 

0.01115015 

0.00019902 

H 

0.04424930 

0.55072200 

0.24769920 

0.02387134 

Theoretical 

MIN 

0.00000002 

0.00108552 

0.00026016 

0.00000010 

K-S 

AVG 

0.00847473 

0.06912110 

0.02318812 

0.00011734 

Statistics 

MAX 

0.02679400 

0.17418700 

0.07323483 

0.00092128 

Empirical 

MIN 

0.00000000 

0.00025278 

0.00002620 

0.00000000 

K-S 

AVG 

0.00712015 

0.02764040 

0.01552219 

0.00003268 

Statistics 

MAX 

0.05078420 

0.14837300 

0.08982747 

0.00079356 

Experiment  Number:  5B 

Samples:  30  samples  of  50  elements  each 

Original  function:  lambda  approximation  of  a  gamma  function 

Method  of  fitment:  conventional  method  of  moments 
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Figure  F.2  Plots  of  50-element  samples  from  an  app 
with  conventional  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.57285600 

1.09224000 

0.76092617 

0.01117079 

Moments  of 

2nd 

0.55903500 

1.59726000 

0.92062887 

0.03486003 

Sample  Data 

3rd 

0.62038900 

1.84820000 

1.14948247 

0.08927224 

4th 

2.76508000 

7.84228000 

4.45436900 

1.66842836 

D 

-0.54812700 

0.11526300 

-0.22124468 

0.03425355 

Lambda 

0.00056224 

0.25148400 

0.13620528 

0.00450865 

Parameters 

0.00000012 

0.04750490 

0.00918425 

0.00014064 

H 

0.00049871 

0.39069200 

0.17601228 

0.00988754 

Theoretical 

MIN 

0.00000137 

0.00071150 

0.00012467 

0.00000004 

K-S 

AVG 

0.00693039 

0.04841000 

0.01740729 

0.00007832 

Statistics 

MAX 

0.01557750 

0.13230700 

0.05959370 

0.00089597 

Empirical 

MIN 

0.00000000 

0.00013362 

0.00002324 

0.00000000 

K-S 

AVG 

0.00691078 

0.02198690 

0.01272067 

0.00001652 

Statistics 

MAX 

0.03848610 

0.10369800 

0.07009377 

0.00039885 

Experiment  Number:  5C 

Samples:  30  samples  of  100  elements  each 

Original  function:  lambda  approximation  of  a  gamma  function 

Method  of  fitment:  conventional  method  of  moments 
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Figure  F.3  Plots  of  100-element  samples  from  an  approximate  gamma  distribution  fitted 
with  conventional  moments 
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Appendix  G.  Experiment  Results:  Exponential  Distribution ,  Conventional 

Moments 

Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.67087600 

1.36818000 

0.93589840 

0.02967292 

Moments  of 

2nd 

0.33428800 

2.25004000 

0.90130390 

0.17730260 

Sample  Data 

3rd 

0.59918000 

2.66973000 

1.37864207 

0.30437496 

4th 

2.14681000 

11.12120000 

4.60744700 

5.42776605 

-0.37838800 

1.85005000 

0.11786311 

0.24590210 

Lambda 

I 

-0.04330980 

0.94384000 

0.15224372 

0.03758558 

Parameters 

-0.00059758 

2.73062000 

0.17564726 

0.44691656 

N 

-0.03724520 

0.74914900 

0.19003573 

0.04033796 

Theoretical 

MIN 

0.00000012 

0.00091408 

0.00035097 

0.00000012 

K-S 

AVG 

0.00890602 

0.05396200 

0.02693338 

0.00013563 

Statistics 

MAX 

0.03926670 

0.19467400 

0.11733060 

0.00155823 

Empirical 

MIN 

0.00000000 

0.00057285 

0.00004606 

0.00000001 

K-S 

AVG 

0.00919480 

0.03766350 

0.01982872 

0.00005756 

Statistics 

MAX 

0.05754620 

0.19870600 

— 

0.13499841 

0.00137565 

. 

Experiment  Number:  6A 

Samples:  30  samples  of  25  elements  each 

Original  function:  lambda  approximation  of  am  exponential  function 
Method  of  fitment:  conventional  method  of  moments 
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Figure  G.l  Plots  of  25-element  samples  from  an  approximate  exponential  distribution  fitted 
with  conventional  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Statistic 


Minimum 

Maximum 

Average 

Variance 

0.75096200 

1.44070000 

0.95519720 

0.01969888 

0.42602100 

2.27778000 

0.86712400 

0.12255268 

0.75023400 

2.02712000 

1.38619813 

E 

.14392488 

2.41646000 

8.46280000 

4.84441167 

2.83071526 

-0.32591000 

0.78509000 

0.13776838 

0.10003427 

-0.00229114 

0.27153700 

0.09106265 

0.01129436 

-0.00066755 

0.00000306 

-0.00014337 

E 

.00000003 

-0.00178665 

0.45000500 

0.13092729 

E 

.02572388 

0.00000578 

0.00186403 

0.00032063 

0.00000018 

0.00646737 

0.07149220 

0.02238879 

0.00013138 

0.04346470 

0.17194900 

0.09520559 

0.00112353 

0.00000000 

0.00019934 

0.00002672 

E 

.00000000 

0.00599659 

0.02537560 

0.01523655 

0.00002695 

0.05326360 

0.18218700 

0.10607955 

0.00109799 

Experiment  Number:  6B 

Samples:  30  samples  of  50  elements  each 

Original  function:  lambda  approximation  of  an  exponential  function 
Method  of  fitment:  conventional  method  of  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


"• .  -  - - - - 

Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.76118800 

1.29035000 

0.94777557 

0.01100065 

Moments  of 

2nd 

0.44768600 

1.77345000 

0.88839270 

0.05455836 

Sample  Data 

3rd 

1.05478000 

2.40918000 

1.59824167 

0.13811616 

4th 

3.27297000 

11.77830000 

5.90627500 

4.32608383 

D 

-0.24337700 

0.62057600 

0.14862391 

0.06535422 

Lambda 

H 

-0.05348860 

0.22742200 

- 

0.05582116 

0.00799021 

Parameters 

-0.00045597 

0.00000376 

-0.00010683 

0.00000002 

H 

-0.04584390 

0.29532000 

0.07217256 

0.01275286 

Theoretical 

MIN 

0.00000022 

0.00091408 

0.00015843 

0.00000006 

K-S 

AVG 

0.00800442 

0.04840770 

0.01663250 

0.00006246 

Statistics 

MAX 

0.02207040 

0.12617800 

0.07962476 

0.00095739 

Empirical 

MIN 

0.00000000 

0.00019789 

0.00002820 

0.00000000 

K-S 

AVG 

0.00508115 

0.02284340 

0.01257003 

0.00001738 

Statistics 

MAX 

0.03472890 

0.14404300 

0.08528730 

0.00082612 

Experiment  Number:  6C 

Samples:  30  samples  of  100  elements  each 

Original  function:  lambda  approximation  of  am  exponential  function 
Method  of  fitment:  conventional  method  of  moments 
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Figure  G.3  Plots  of  100-element  samples  from  an  approximate  exponential  distribution 
fitted  with  conventional  moments 
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Appendix  H.  Experiment  Results:  Normal  Distribution,  Alternate  Moments 

Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

-0.32920000 

0.33291600 

-0.08192267 

0.03075791 

Moments  of 

2nd 

0.68745300 

1.54415000 

1.03898053 

0.05250603 

Sample  Data 

3rd 

0.46901900 

1.74990000 

1.01318137 

0.12075743 

4th 

1.87260000 

3.05616000 

2.42375433 

0.07598671 

D 

-1.81571000 

2.08373000 

0.01490352 

0.78299881 

Lambda 

-0.15520400 

0.40728500 

0.21428818 

0.01661824 

Parameters 

-0.06558470 

0.84846300 

0.22445681 

0.04732022 

N 

-0.08777880 

0.61003900 

0.18532574 

0.04073945 

Theoretical 

MIN 

0.00000194 

0.00132908 

0.00014137 

0.00000008 

K-S 

AVG 

0.01242070 

0.05460580 

0.03524031 

0.00015237 

Statistics 

MAX 

0.03273380 

0.15955000 

0.09805817 

0.00119411 

Empirical 

MIN 

0.00000000 

0.00013531 

0.00001919 

0.00000000 

K-S 

AVG 

0.01167520 

0.02929620 

0.02041950 

0.00002269 

Statistics 

MAX 

0.05895210 

0.14434900 

0.09787561 

0.00047475 

Experiment  Number:  7  A 

Samples:  30  samples  of  25  elements  each 

Original  function:  lambda  approximation  of  a  normal  function 

Method  of  fitment:  alternate  method  (Q3,  Q4) 
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Figure  H.l  Plots  of  25-element  samples  from  an  approximate  normal  distribution  fitted 
with  alternate  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

-0.24640400 

0.34302800 

-0.04351340 

0.01692519 

Moments  of 

2nd 

0.67163200 

1.32877000 

1.01043343 

0.02903252 

Sample  Data 

3rd 

0.61927500 

1.38304000 

0.97787483 

0.03994153 

4th 

1.98703000 

2.95123000 

2.44442767 

0.04833684 

D 

-1.70150000 

1.52504000 

-0.03031136 

0.54837742 

Lambda 

-0.08009890 

0.35489200 

0.22881436 

0.00785004 

Parameters 

-0.03850840 

0.70451800 

0.19743541 

0.02759746 

H 

-0.03574680 

0.76673800 

0.20614792 

0.04422529 

Theoretical 

MIN 

0.00000051 

0.00128625 

0.00018631 

0.00000014 

K-S 

AVG 

0.01011840 

0.05588620 

0.02453645 

0.00011004 

Statistics 

MAX 

0.02602540 

0.12451500 

0.06674605 

0.00087331 

Empirical 

MIN 

0.00000000 

0.00009787 

0.00002051 

0.00000000 

K-S 

AVG 

0.00994767 

0.02185460 

0.01436887 

0.00001010 

Statistics 

MAX 

0.04597370 

0.13449200 

0.07069398 

0.00036417 

Experiment  Number:  7B 

Samples:  30  samples  of  50  elements  each 

Original  function:  lambda  approximation  of  a  normal  function 

Method  of  fitment:  alternate  method  ( Q3 ,  Q±) 
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with  alternate  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

-0.21766800 

0.23810900 

-0.04653627 

0.01058862 

Moments  of 

2nd 

0.79556500 

1.19861000 

0.99209797 

0.00903663 

Sample  Data 

3rd 

0.65539800 

1.23400000 

0.98769230 

0.02105896 

4th 

2.13103000 

2.89916000 

2.50212733 

0.02084603 

n 

-1.15944000 

0.85998900 

-0.04724554 

0.17590329 

Lambda 

-0.02939990 

0.34285700 

0.23702967 

0.00653531 

Parameters 

-0.01272650 

0.38912000 

0.18013660 

0.00888541 

H 

-0.01722070 

0.69886500 

0.19117286 

0.01827860 

Theoretical 

MIN 

0.00000028 

0.00127119 

0.00015771 

0.00000012 

K-S 

AVG 

0.00778402 

0.03931520 

0.01862765 

0.00007079 

Statistics 

MAX 

0.01747080 

0.09863610 

0.05115571 

0.00055662 

Empirical 

MIN 

0.00000000 

0.00005997 

0.00001558 

0.00000000 

K-S 

AVG 

0.00640562 

0.02018290 

0.01038700 

0.00000838 

Statistics 

MAX 

0.03170860 

0.08823640 

0.05005300 

0.00024256 

Experiment  Number:  7C 

Samples:  30  samples  of  100  elements  each 

Original  function:  lambda  approximation  of  a  normal  function 

Method  of  fitment:  alternate  method  (Q3,  Q4) 
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Appendix  I.  Experiment  Results:  Gamma  Distribution ,  Alternate  Moments 

Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.48061600 

1.17188000 

0.74208173 

0.03016900 

Moments  of 

2nd 

0.44652400 

1.83663000 

0.93963863 

0.11671945 

Sample  Data 

3rd 

1.13009000 

4.67810000 

2.52182233 

0.88282300 

4th 

1.85934000 

3.51239000 

2.55625267 

0.13492001 

D 

-0.71157300 

0.41072000 

-0.20031334 

0.07819155 

Lambda 

1 

-0.32067000 

0.42483700 

0.13612248 

0.02632154 

Parameters 

-0.03314990 

0.14390200 

0.01575699 

0.00103320 

K9 

-0.17998800 

0.57410200 

0.20646465 

0.04281049 

Theoretical 

MIN 

0.00000011 

0.00088285 

0.00015506 

0.00000006 

K-S 

AVG 

0.01037650 

0.05087940 

0.02823904 

0.00010827 

Statistics 

MAX 

0.02847320 

0.16217900 

0.09367535 

0.00099154 

Empirical 

MIN 

0.00000000 

0.00049341 

0.00003267 

0.00000001 

K-S 

AVG 

0.01000580 

0.03048770 

0.01731186 

0.00002533 

Statistics 

MAX 

0.06234880 

0.16882500 

0.10338170 

0.00073998 

Experiment  Number:  8A 

Samples:  30  samples  of  25  elements  each 

Original  function:  lambda  approximation  of  a  gamma  function 
Method  of  fitment:  alternate  method  (Q 3,  Q 4) 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 


1st 

Moments  of  2nd 
Sample  Data  3rd 


Lambda 

Parameters 


Theoretical 

K-S 

Statistics 


Empirical 

K-S 

Statistics 


Statistic 

Minimum 

Maximum 

Average 

Variance 

0.55609200 

1.23099000 

0.76728553 

0.01912877 

0.51725800 

1.95520000 

0.90982420 

0.08046629 

1.53465000 

3.49965000 

2.45939067 

0.31804869 

2.09833000 

3.08585000 

2.54661433 

0.07307762 

-0.59665300 

0.20914800 

-0.16800964 

0.05185907 

-0.18378400 

0.34113400 

0.13081282 

0.01268965 

-0.03185830 

0.04747990 

0.00908726 

0.00021125 

-0.10452200 

0.51679700 

0.18775621 

0.02943946 

0.00000032 

0.00088285 

0.00018192 

0.00000008 

0.00824128 

0.05872040 

0.02087910 

0.00009499 

0.02166040 

0.12191200 

0.06534404 

0.00081286 

0.00000000 

0.00006081 

0.00000854 

_ 

0.00000000 

0.00676882 

0.01756250 

0.01199648 

0.00000938 

0.04512980 

0.12361500 

0.06964856 

0.00033589 

Experiment  Number:  8B 

Samples:  30  samples  of  50  elements  each 

Original  function:  lambda  approximation  of  a  gamma  function 

Method  of  fitment:  alternate  method  (Q3)  Qt) 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.57285600 

1.09224000 

0.76092617 

0.01117079 

Moments  of 

2nd 

0.55903500 

1.59726000 

0.92062887 

0.03486003 

Sample  Data 

3rd 

1.62456000 

3.39448000 

2.53449067 

0.18600054 

4th 

2.24375000 

3.27435000 

2.66284800 

0.04042513 

-0.50730100 

0.16275400 

-0.09489049 

0.02077248 

Lambda 

-0.29228700 

0.22253200 

0.09354151 

0.01083183 

Parameters 

-0.03866660 

0.03801980 

0.00958450 

0.00019587 

H 

-0.17385200 

0.33751500 

0.11798481 

0.01265736 

Theoretical 

MIN 

0.00000053 

0.00088285 

0.00011864 

0.00000006 

K-S 

AVG 

0.00682230 

0.04039520 

0.01554683 

0.00006200 

Statistics 

MAX 

0.01639940 

0.09027510 

0.04910512 

0.00045255 

Empirical 

MIN 

0.00000000 

0.00005496 

0.00001453 

0.00000000 

K-S 

AVG 

0.00507575 

0.01403150 

0.00871258 

0.00000485 

Statistics 

MAX 

0.03017440 

0.08858540 

0.04963251 

0.00022438 

Experiment  Number:  8C 

Samples:  30  samples  of  100  elements  each 

Original  function:  lambda  approximation  of  a  gamma  function 

Method  of  fitment:  alternate  method  ( Q3 ,  Q 4) 
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Figure  1.3  Plots  of  100-element  samples  from  an  approximate  gamma  distribution  fitted 
with  alternate  moments 
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Appendix  J.  Experiment  Results:  Exponential  Distribution,  Alternate 

Moments 

Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.67087600 

1.36818000 

0.93589840 

0.02967292 

Moments  of 

2nd 

0.33428800 

2.25004000 

0.90130390 

0.17730260 

Sample  Data 

3rd 

2.21996000 

7.41092000 

4.21758600 

2.46529077 

4th 

1.85761000 

3.74385000 

2.60321133 

0.20674032 

D 

-0.22480400 

0.57029400 

-0.01852245 

0.02455917 

Lambda 

-0.49614100 

0.35312500 

0.05022049 

0.03538417 

Parameters 

-0.04763650 

0.00964982 

-0.00115794 

0.00008390 

H 

-0.30838200 

0.40841600 

0.09051723 

0.03646706 

Theoretical 

MIN 

0.00000046 

0.00872040 

0.00057302 

0.00000248 

K-S 

AVG 

0.00395157 

0.04919680 

0.02529038 

0.00012810 

Statistics 

MAX 

0.03431670 

0.23971400 

0.08873605 

0.00173606 

Empirical 

MIN 

0.00000000 

0.00029412 

0.00004913 

0.00000001 

K-S 

AVG 

0.00947864 

0.02971550 

0.01766614 

0.00003414 

Statistics 

MAX 

0.06381820 

0.23657200 

0.11254081 

0.00130056 

Experiment  Number:  9 A 

Samples:  30  samples  of  25  elements  each 

Original  function:  lambda  approximation  of  an  exponential  function 
Method  of  fitment:  alternate  method  (Qa,  Qi) 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.75096200 

1.44070000 

0.95519720 

0.01969888 

Moments  of 

2nd 

0.42602100 

2.27778000 

0.86712400 

0.12255268 

Sample  Data 

3rd 

2.54713000 

6.81954000 

4.17582967 

1.11603068 

4th 

2.09800000 

3.12824000 

2.58020400 

0.09538596 

-0.15976300 

0.20232000 

-0.00059699 

0.00590568 

Lambda 

-0.18197200 

0.27843100 

0.05597127 

0.01261634 

Parameters 

1 

-0.00878814 

0.00441249 

0.00000719 

0.00000366 

H 

-0.11115800 

0.33144400 

0.07108713 

0.01512170 

Theoretical 

MIN 

0.00000022 

0.00183874 

0.00040423 

0.00000019 

K-S 

AVG 

0.00570432 

0.06348880 

0.01867957 

0.00012561 

Statistics 

MAX 

0.02510570 

0.13498200 

0.05966165 

0.00081784 

Empirical 

MIN 

0.00000000 

0.00014532 

0.00003065 

0.00000000 

K-S 

AVG 

0.00585411 

0.02032030 

0.01239072 

0.00001365 

Statistics 

MAX 

0.04612920 

0.11354700 

0.07479089 

0.00035223 

Experiment  Number:  9B 

Samples:  30  samples  of  50  elements  each 

Original  function:  lambda  approximation  of  an  exponential  function 
Method  of  fitment:  alternate  method  (Q3,  Q4) 
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Figure  J.2  Plots  of  50-element  samples  from  an  approximate  exponential  distribution  fitted 
with  alternate  moments 
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Monte  Carlo  Experiment  Results 
Statistical  Summary 


Experiment 

Result 

Statistic 

Minimum 

Maximum 

Average 

Variance 

1st 

0.76118800 

1.29035000 

0.94777557 

0.01100065 

Moments  of 

2nd 

0.44768600 

1.77345000 

0.88839270 

0.05455836 

Sample  Data 

3rd 

3.01858000 

6.07224000 

4.27603400 

0.65647523 

4th 

2.29008000 

3.43683000 

2.72610767 

0.05867238 

Ai 

-0.09161470 

0.11692500 

-0.01957376 

0.00235948 

Lambda 

Aa 

-0.09522080 

0.17726900 

0.04257708 

0.00597940 

Parameters 

As 

-0.00279383 

0.00386742 

0.00007093 

0.00000107 

A4 

-0.08411940 

0.22178400 

0.05007467 

0.00660096 

Theoretical 

MIN 

0.00000203 

0.00091407 

0.00026584 

0.00000009 

K-S 

AVG 

0.00306969 

0.04246320 

0.01428211 

0.00006877 

Statistics 

MAX 

0.01800200 

0.09584070 

0.04437378 

0.00042776 

Empirical 

MIN 

0.00000000 

0.00011963 

0.00002374 

0.00000000 

K-S 

AVG 

0.00452803 

0.01687420 

0.00906632 

0.00000789 

Statistics 

MAX 

0.02731380 

0.10242200 

0.05183214 

0.00040186 

Experiment  Number:  9C 

Samples:  30  samples  of  100  elements  each 

Original  function:  lambda  approximation  of  an  exponential  function 
Method  of  fitment:  alternate  method  ( Q3 ,  Q 4) 
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Figure  J.3  Plots  of  100-element  samples  from  an  approximate  exponential  distribution  fit¬ 
ted  with  alternate  moments 
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Appendix  K.  FORTRAN  Programs:  Random  Variate  Data  Sample 

Generator 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMEHTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  SAMPGEN 

*  SUBROUTINE  PURPOSE:  GENERATE  A  FILE  OF  RANDOM  VARIATES  FROM  THE  GLD 

*  PERCENTILE  FUNCTION 

*  '  ROUTINE  DESCRIPTION: 

*  This  program  loads  an  array  with  a  random  number  generator,  generating 

*  pseudo-random  numbers  from  zero  to  one.  Each  random  variate  is  then 

*  converted  to  a  GLD  function  using  the  GLD  quantile  function,  with 

*  the  lambda  values  set  as  parameters  of  the  program. 

*  The  parameter  section  of  the  program  provides  the  programmer  with 

*  several  sets  of  distribution  parameters  to  choose  from.  The  programner 

*  chooses  the  desired  distribution  by  commenting  out  all  the  other 

*  parameter  statements.  For  example ,  in  this  copy,  the  NORMAL  distribution 

*  is  selected. 

*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  OCTOBER  1993 

*  INPUT  PARAMETERS:  HOME. 

*  OUTPUT  PARAMETERS:  NONE. 

******* *************************** *♦♦*♦*♦* *************************** ******** 
SUBROUTINE  SAMPGEN 

******  DECLARE  VARIABLES,  PARAMETERS,  AND  INTRINSIC  FUNCTIONS 
DOUBLE  PRECISION  CHART  (  1:1000  ),  SAMPLE  (  1:1000  ) 

DOUBLE  PRECISION  LI,  L2,  L3,  L4,  P 
CHARACTER  ANSWER 
INTEGER  HUM,  I 

******  NORMAL  DISTRIBUTION  PARAMETERS 

******  mean=0 ,  sigma  2=1,  skewness=0,  kurtosis=3 

PARAMETER  (  L1=0.0D0,  L2=0.19750D0,  L3=0.13490D0,  L4=0.13490D0  ) 

******  UNIFORM  DISTRIBUTION  PARAMETERS 

******  mean=0 ,  sigma  2=1,  skewness=0,  kurtosis=1.8 

**  PARAMETER  (  L1=0.0D0,  L2=0 . 57740D0 ,  L3=1.0D0,  L4=1.0D0  ) 
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******  EXPONENTIAL  DISTRIBUTION  PARAMETERS 
*****41  mean=0 ,  sigma  2=1,  skewness=2,  kurtosis=9 

**  PARAMETER  (Ll=-0.9930D0,L2=-0. 001081, L3=-0. 0000041,  L4=-0. 00108) 

******  STUDENT'S  T  DISTRIBUTION  PARAMETERS 

******  mean=0 ,  sigma  2=1,  skevness=0,  kurtosis=9 

**  PARAMETER  (L1=0.0D0,L2=-0.32030D0,L3=-0.13590D0,L4=-0.13590D0) 

******  VEIBULL  DISTRIBUTION  PARAMETERS 

******  mean=0,  sigma  2=1,  ske«ness=l .4,  kurtosis=5.8 

**  PARAMETER  (L1=-0.8440DO,L2=0.05380DO,L3=0.006530DO,L4=0.05470DO) 

******  GAMMA  DISTRIBUTION  PARAMETERS 

******  mean=0,  sigma  2=1,  skewness=l,  kurtosis=4.6 

**  PARAMETER  (L1=-0.6380D0,L2=0.07410D0,L3=0.01820D0,L4=0.06970D0) 

******  beta  DISTRIBUTION  PARAMETERS 

******  mean=0,  sigma  2=1,  skevness=l .7,  kurtosis=5.2 

**  PARAMETER  (  L1=0.0D0,  L2=0.0D0,  L3=0.0D0,  L4=O.ODO  ) 

INTRINSIC  MOD 

******  PROMPT  USER  FOR  INPUTS 
300  PRINT  *,  '  ' 

PRINT  *,  'DO  YOU  WANT  TO  GENERATE  A  DATA  SAMPLE  FILE?' 

READ  *,  ANSWER 

IF  ((ANSWER. NE.'Y'). AND. (ANSWER. NE.'y'))  THEN 
GOTO  200 
END  IF 

PRINT  *,  '  ' 

PRINT  *,  'ENTER  THE  NUMBER  OF  ELEMENTS  IN  THE  SAMPLE:' 

READ  *,  NUM 
PRINT  *,  '  ' 

PRINT  *,  'ENTER  THE  SEED  FOR  THE  RANDOM  NUMBER  GENERATOR: ' 

PRINT  *,  '  ' 

PRINT  *,  '(AN  INTEGER  BETWEEN  2  AND  2147483647)' 

READ  *,  SEED 
PRINT  *,  '  ' 

******  GENERATE  ARRAY  OF  RANDOM  NUMBERS 
CALL  DURAND  (SEED,  NUM,  CHART) 

******  GENERATE  ARRAY  OF  GLD  RANDOM  VARIATES 
DO  100  I  =  1,  NUM 
P  =  CHART  (I) 

SAMPLE  (I)  =  LI  +  ( (P**L3) - ( (1-P) **L4) ) /L2 
100  CONTINUE 

CALL  SAMPFILE  (NUM,  SAMPLE) 

GOTO  300 
200  RETURN 
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END 


***************************************************************************** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  DURAND 

*  SUBROUTINE  PURPOSE:  PSEUDO-RANDOM  NUMBER  GENERATOR 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  uses  a  multiplicative  congruential  generator  with  a 

*  base  of  2**31-1  and  a  multiplier  of  7**5. 

*  AUTHOR(S) : 

*  J.  R.  M.  HOSKING 

*  IBM  RESEARCH  DIVISION 

*  T.  J.  VATSON  RESEARCH  CENTER 

*  YORKTOHN  HEIGHTS 

*  NEV  YORK  10598,  USA 

*  DATE:  JULY  1988 

*  INPUT  PARAMETERS: 

*  SEED  *  IN/OUT*  SEED  FOR  RANDOM  NUMBER  GENERATOR.  SHOULD  BE  A  WHOLE 

*  NUMBER  IN  THE  RANGE  2D0  TO  2147483647D0. 

*  N  *  INPUT*  NUMBER  OF  NUMBERS  TO  BE  GENERATED 

*  OUTPUT  PARAMETERS: 

*  I  *OUTPUT*  ARRAY  OF  LENGTH  N.  ON  EXIT,  CONTAINS  RANDOM  NUMBERS. 

***************************************************************************** 

SUBROUTINE  DURAND ( SEED, N,X) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  X(N) 

DATA  AMULT , BASE/16807D0 , 2147483647D0/ 

DO  10  1=1, N 

SEED=DMOD (SEED* AMULT , BASE) 

X(I)=SEED/BASE 
10  CONTINUE 
RETURN 
END 

it**************************************************************************** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE : SAMPFILE 

*  SUBROUTINE  PURPOSE:  OUTPUT  ARRAY  OF  RANDOM  SAMPLE  DATA  ELEMENTS  TO  FILE 
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*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  inputs  an  array  of  random  variates  from  the  GLD  random 

*  8 ample  generator  and  loads  them  into  an  output  file.  The  file  name  is 

*  input  by  the  user  when  prompted  by  this  routine. 

*  AUTHOR(S) : 

*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  OCTOBER  1993 

*  INPUT  PARAMETERS: 

*  HUM:  INTEGER  PARAMETER  number  of  points  in  set 

*  LAMBDA:  DOUBLE  PRECISION  PARAMETER  data  array 

*  NAME:  CHARACTER  STRING  output  file  name 

*  OUTPUT  PARAMETERS:  NONE. 

***************************************************************************** 
SUBROUTINE  SAMPFILE  (  BUM,  LAMBDA  ) 

******  DECLARE  LOCAL  VARIABLES 

DOUBLE  PRECISION  LAMBDA  (  1:1000  ) 

INTEGER  BUM,  I,  II 
CHARACTER* 12  NAME 

******  PROMPT  USER  FOR  NAME  OF  FILE 

PRINT  *,  ’ENTER  THE  NAME  OF  THE  OUTPUT  FILE: ’ 

READ  (*,  ’ (A12) ’  )  NAME 

******  open  OUTPUT  FILE  AND  DOWNLOAD  FORMATTED  DATA 

OPEN  (UNIT=30 ,  FILE=NAME,  STATUS® ’NEW’ ,  IOSTAT=II,  ERR=200) 

10  FORMAT  (F9.4) 

******  download  array 
DO  100  I  -  1,  BUM 

WRITE  (30,  10)  LAMBDA  (I) 

100  CONTINUE 

******  close  output  file 
CLOSE  (30) 

******  RETURN  TO  MAIN  PROGRAM 
RETURN 

******  ERROR  TRAP 

200  PRINT  *,  ’CANNOT  OPEN  FILE  ’,  NAME,  ’,  ERROR®  ’,  II 
RETURN 
END 
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Appendix  L.  FORTRAN  Programs:  Data  Sample  Input  and  Sort  Routines 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IH  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L -MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  POINTSET 

*  SUBROUTINE  PURPOSE:  Input  a  sample  data  file  into  an  array 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  prompt s  the  user  for  the  name  of  the  input  file,  opens  the 

*  data  input  file,  and  retreives  the  set  of  points  that  make  up  the  sample 

*  data.  It  then  loads  the  points  into  an  array.  The  array  is  limited  to 

*  1000  points. 

*  AUTHOR(S) : 

*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  HRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  OCTOBER  1993 

*  INPUT  PARAMETERS:  HONE. 

*  OUTPUT  PARAMETERS: 

*  CHART  -  1000  element  array  containing  points  from  sample  data  file. 

*  NUM  -  integer  containing  number  of  elements  in  the  CHART  array. 
***************************************************************************** 

SUBROUTINE  POINTSET  (  CHART,  NUM  ) 

******  DECLARE  LOCAL  VARIABLES 
INTEGER  NUM,  I,  II 
DOUBLE  PRECISION  CHART  (  1 : 1000  ) 

CHARACTER* 12  NAME 

******  PROMPT  USER  FOR  NAME  OF  FILE 
PRINT  *,  ’  » 

PRINT  *,  ’ENTER  THE  NAME  OF  THE  SAMPLE  DATA  FILE:’ 

READ  (*,  ’ (A12) ’  )  NAME 
******  OPEN  INPUT  FILE  AND  RETREIVE  DATA 

OPEN  (UNIT=10 ,  FILE=NAME,  STATUS=’ OLD’ ,  IOSTAT=II,  ERR=130) 

******  READ  ELEMENTS  INTO  ARRAY 
NUM  =  0 

DO  110  I  =  1,  1000 

READ  (10,100,  IOSTAT=II,  ERR=140,  END=120)  CHART  (I) 

100  FORMAT  (F9.4) 

NUM  =  NUM  +  1 
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110  COBTIHUE 

PRIHT  *,  '1000  ELEMENTS  READ  INTO  ARRAY' 

******  CLOSE  INPUT  FILE 
120  CLOSE  (10) 

******  RETURN  TO  MAIN  PROGRAM 
RETURN 

******  ERROR  TRAP 

130  PRINT  *,  'CANNOT  OPEN  FILE  ',  NAME,  ',  ERROR*  ',  II 
RETURN 

140  PRINT  *,  'CANNOT  READ  FILE  ',  NAME,  ',  ERROR*  ',  II 
RETURN 
END 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER'S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  SORT 

*  SUBROUTINE  PURPOSE:  SORT  AN  ARRAY  INTO  ASCENDING  ORDER 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  sorts  the  input  array  of  real  numbers  into  ascending  order 

*  and  outputs  the  sorted  err  ay .  It  does  not  use  a  temporary  array.  The 

*  method  used  is  a  shell  sort  with  a  sequence  of  increments. 

*  AUTHOR(S) : 

*  J.  R.  M.  HOSKING 

*  IBM  RESEARCH  DIVISION 

*  T.  J.  MATSON  RESEARCH  CENTER 

*  YORKTOVN  HEIGHTS 

*  NEW  YORK  10598,  USA 

*  DATE:  JULY  1988 

*  INPUT  PARAMETERS: 

*  X  -  Array  of  N  elements,  double  precision  real  numbers,  contains  the 

*  sample  data  elements  to  be  sorted. 

*  N  -  Integer  containing  the  number  of  elements  in  the  data  sample  array. 

*  OUTPUT  PARAMETERS : 

*  X  -  Array  of  N  elements,  double  precision  real  numbers,  contains  the 

*  sorted  sample  data  elements. 

*  N  -  Integer  containing  the  number  of  elements  in  the  data  sample  array. 
***************************************************************************** 

SUBROUTINE  SORT  (X,N) 

******  DECLARE  VARIABLES 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
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DOUBLE  PRECISION  1(H) 

******  ERROR  CHECKIHG 
IF (H.LE.l) RETURN 
J-4 

DO  10  1=1,100 
J«3*J+1 

IF(J.GE.H)GOTO  20 
10  CONTINUE 
20  CONTINUE 
M»(J/3) 

DO  60  MM=1,100 
M=M/3 

IF (M.EQ.O) RETURN 
DO  60  I=M+1,N 
TEST=X(I) 

J=I 

DO  30  JJ=1,100 
J=J-M 

IF(J.LE.O)GOTO  40 
IF (TEST. GE.X(J)) GOTO  40 
30  X(J+M)*X(J) 

40  CONTINUE 

60  X(J+M)=TEST 
60  CONTINUE 
END 
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Appendix  M.  FORTRAN  Program:  Sample  L-moment  Computer 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA"* 

*  SUBROUTINE  TITLE:  SAMLMR 

*  SUBROUTINE  PURPOSE:  COMPUTE  THE  SAMPLE  L-MOMENTS  OF  A  SAMPLE  DATA  ARRAY 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  receives  the  number  of  L-moment s  desired  (up  to  20)  ft  array 

*  of  sample  data  and  computes  the  L-moment s  for  that  sample  data.  It  begins 

*  by  computing  the  order  statistics  and  applying  the  plotting  positions. 

*  Then  it  computes  the  L-moment  functions  of  the  order  statistics . 

*  FOR  UNBIASED  ESTIMATES  (OF  THE  LAMBDA’S)  SET  A=B=ZER0 .  OTHERWISE, 

*  PLOTTING-POSITION  ESTIMATORS  ARE  USED,  BASED  ON  THE  PLOTTING  POSITION 

*  ( J+A) / (N+B)  FOR  THE  J’TH  SMALLEST  OF  N  OBSERVATIONS.  FOR  EXAMPLE, 

*  A--0.35D0  AND  B=0 . ODO  YIELDS  THE  ESTIMATORS  RECOMMENDED  BY 

*  HOSKING  ET  AL.  (1985,  TECHNOMETRICS)  FOR  THE  GEV  DISTRIBUTION. 

*  AUTHOR(S) : 

*  J.  R.  M.  HOSKING 

*  IBM  RESEARCH  DIVISION 

*  T.  J.  WATSON  RESEARCH  CENTER 

*  YORKTOWN  HEIGHTS 

*  NEW  YORK  10598,  USA 

*  DATE:  JULY  1988 

*  INPUT  PARAMETERS: 

*  X  *  INPUT*  ARRAY  OF  LENGTH  N.  CONTAINS  THE  DATA,  IN  ASCENDING  ORDER. 

*  N  *  INPUT*  NUMBER  OF  DATA  VALUES 

*  NMOM  *  INPUT*  NUMBER  OF  L-MOMENTS  TO  BE  FOUND.  AT  MOST  MAX (N, 20)  . 

*  A  *  INPUT*  )  PARAMETERS  OF  PLOTTING 

*  B  *  INPUT*  )  POSITION  (SEE  ABOVE) 

*  OUTPUT  PARAMETERS: 

*  XMOM  *OUTPUT*  ARRAY  OF  LENGTH  NMOM.  ON  EXIT,  CONTAINS  THE  SAMPLE 

*  L-MOMENTS  L-l,  L-2,  T-3,  T-4,  ...  . 
***************************************************************************** 

SUBROUTINE  SAMLMR (X,N, XMOM, NMOM, A, B) 

******  DECLARE  VARIABLES  AND  INITIALIZE  PARAMETERS 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  X(N) ,XMOM(NMOM) ,SUM(20) 

DATA  ZERO /ODO/ , ONE/ IDO/ 

******  ERROR  CHECKING:  Too  many  L-moments  for  routine  or  sample  size 
IF (NMOM . GT . 20 . OR . NMOM . GT . N) GOTO  1000 
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******  IHITIALIZE  ARRAY  VALUES  TO  ZERO 
DO  10  J-l.NMOM 
10  SUM (J) -ZERO 

******  BYPASS  PLOTTING-POSITION  ROUTIME  IF  A  AND  B  ARE  ZERO 
IF (A. EQ. ZERO. AND. B.EQ. ZERO) GOTO  50 
******  ERROR  CHECKING:  A  <=  -1  or  a  >a  B 
IF(A.LE.-ONE.OR. A.GE.B)GOTO  1010 
*  PLOTTING-POSITION  ESTIMATES  OF  PWM’S 


* 


*  PLOTTING-POSITION  (PPOS) 

* 

* 

* 

DO  30  1*1, N 

* 

PPOS**  ( I+A)  /  (N+B) 

* 

TERM=X ( I ) 

SUM  ( 1 )  =SUM  ( 1 )  +TERM 

* 

DO  20  J**2,NM0M 

* 

TERM=TERM*PPOS 

* 

20  SUM(J)=SUM(J)+TERM 

* 

* 


Ith  Value  of  N  values  + 

Position  Adjustment  A  / 

N  (sample  size)  + 

Position  Adjustment  B 

for  each  element  of  the  sample, 

compute  the  plotting  position, 

compute  the  sum  for  the  mean  (L-l)  . 

for  second  L -moment  and  higher, 

the  element  z  plotting  position, 

add  the  term  to  the  running  total 
for  the  Jth  L -moment . 


30  CONTINUE 

******  DIVIDE  EACH  SUM  VALUE  BY  NUMBER  OF  VALUES  IN  ARRAY  (AVERAGE) 
DO  40  J=1 ,NM0M 
40  SUM(J)=SUM(J)/N 
******  SKIP  NEXT  SECTION 
GOTO  100 

******  UNBIASED  ESTIMATES  OF  PWM'S  (NOT  USED  WITH  GLD) 

50  DO  70  1*1, N 
Z=I 


TERM=X(I) 

SUM  ( 1 )  =SUM  ( 1 )  +TERM 
DO  60  J=2,NM0M 
Z=Z-ONE 


TERM=TERM*Z 

60  SUM( J)=SUM( J)+TERM 
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70  comiuE 

Y=N 

Z-i 

SUM(1)=SUM(1)/Z 
DO  80  J-2.NM0M 
Y-Y-ONE 


Z=Z*Y 


80  SUM(J)=SUM(J)/Z 

******  L-MOMENTS:  Compute  functions  of  PWM's 
100  K-NMOM 
PO-OHE 

IF(NMOM-NMOM/2*2 .EQ .  1)PO=-OHE 
DO  120  KK=2 ,HMOM 

*  for  second  L -moment  and  higher, 

AK-K 
PO=-PO 


* 

* 


* 

110 


* 

* 


120 


P»PO 

TEMP=P*SUM(1) 
DO  110  1=1, K-l 


AI=I 


load  mean  with  adjusted  sign, 
for  remaining  L-moments, 


P=-P* (AK+AI-OHE) * (AK-AI) / (AI*AI) 

compute  function  coefficients, 

TEMP=TEMP+P*SUM( 1+1 ) 

add  to  running  total  of  linear 
function. 

SUM(K)=TEMP 


K=K-1 


******  COMPUTE  MEAN  OF  SAMPLE 
XM0M(1)=SUM(1) 

******  ip  OHLY  1ST  MOMENT  IS  WANTED,  RETURN 
IF (NM0M.EQ.1) RETURN 
XM0M(2)=SUM(2) 

******  ERROR  CHECK  FOR  NON-VARIABLE  SAMPLE  (NO  VARIANCE) 
IF(SUM(2).EQ.ZER0)G0T0  1020 
IF  (NMQM .  EQ .  2)  RETURN 

******  COMPUTE  TAU  VALUES  (TAU-K  =  LAMBDA-K  /  LAMBDA-2) 
DO  130  K=3,NM0M 
130  XMOM(K) =SUM(K) /SUM(2) 

RETURN 

******  ERROR  traps 
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1000  WRITE (6, 7000) 

RETURI 

1010  WRITE (6, 7010) 

RETURI 

1020  WRITE(6,7020) 

RETURI 

7000  FORMAT  (*  ***  ERROR  ***  ROUT  I  ME  SAMLMR  :  PARAMETER  IMOM  IIVALID’) 
7010  FORMAT (’  ***  ERROR  ***  ROUTIIE  SAMLMR 

*  »  PLOTTIIG-POSITIOI  PARAMETERS  IMVALID’) 

7020  FORMAT (»  ***  ERROR  ***  ROUTIIE  SAMLMR  :  ALL  DATA  VALUES  EQUAL’) 
EID 


Appendix  N.  FORTRAN  Program :  L-Moment  to  Lambda  Parameter 

Computer 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IH  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  PELTU4 

*  SUBROUTINE  PURPOSE:  PARAMETER  ESTIMATION  VIA  L-MOMENTS  FOR  THE  4-PARAMETER 

*  GENERALIZED  LAMBDA  DISTRIBUTION 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  receives  the  first  four  sample  L -moments  from  a  data  sample 

*  and  computes  the  four  lambda  parameters  of  the  GLD  using  a  Hevton-Raphson 

*  search  routine  for  solving  two  equations  with  two  unknowns. 

*  IN  GENERAL,  TWO  SETS  OF  PARAMETERS  CAN  BE  FOUND  WHICH  YIELD  A  GIVEN 

*  SET  OF  L-MOMEHTS.  THE  ROUTINE  LOOKS  ONLY  FOR  PARAMETERS  HHICH  SATISFY 

*  -l.LT.BETA.LE.+l  AND  -1 .LT. DELTA. LE.+l .  SUCH  PARAMETERS  DO  HOT  EXIST 

*  FOR  ALL  SETS  OF  L-MOMENTS,  BUT  ARE  UNIQUE  HHEN  THEY  DO  EXIST. 

*  THE  BETA  AND  DELTA  PARAMETERS  ARE  ESTIMATED  USING  NEWTOH-RAPHSON 

*  ITERATION  ON  THE  RELATIONSHIP  BETWEEN  (TAU-3.TAU-4)  AND  (BETA, DELTA) . 

*  THE  CONVERGENCE  CRITERION  IS  THAT  TAU-3  AND  TAU-4  CALCULATED  FROM  THE 

*  ESTIMATED  VALUES  OF  BETA  AND  DELTA  SHOULD  DIFFER  BY  LESS  THAN  IE-6 

*  FROM  THE  VALUES  SUPPLIED  IN  ARRAY  XMOM.  UNLESS  A  PARAMETER  IS  CLOSE 

*  TO  +1 ,  THE  RETURNED  PARAMETER  ESTIMATES  HILL  ALSO  BE  CORRECT  TO 

*  WITHIN  IE-6. 

*  AUTHOR(S) : 

*  J.  R.  M.  HOSKING 

*  IBM  RESEARCH  DIVISION 

*  T.  J.  HATSON  RESEARCH  CENTER 

*  YORKTOWN  HEIGHTS 

*  NEH  YORK  10598,  USA 

*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  OCTOBER  1993 

*  INPUT  PARAMETERS : 

*  XMOM  *  INPUT*  ARRAY  OF  LENGTH  4.  CONTAINS  THE  L-MOMENTS  IN  THE  ORDER 

*  LAMBDA-1,  LAMBDA-2,  TAU-3,  TAU-4. 

*  OUTPUT  PARAMETERS: 


N-l 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


PARA  * OUTPUT*  ARRAY  OF  LEMGTH  5.  CONTAINS  THE  ESTIMATED  PARAMETERS 

II  THE  ORDER  XI,  ALPHA,  BETA,  GAMMA,  DELTA. 

ALPHA  AID  GAMMA  ARE  EQUAL . 

GLD  * OUTPUT*  ARRAY  OF  LEMGTH  4.  COHTAIHS  THE  ESTIMATED  PARAMETERS 

IH  THE  ORDER  LAMBDA-1,  -2,  -3,  -4. 

(THESE  ARE  RAMBERG’S  GLD  PARAMETERS) 

IFAIL  *OUTPUT*  FAIL  FLAG.  OH  EXIT,  IT  IS  SET  TO: 

0  OH  SUCCESSFUL  EXIT 

1  IF  L-MOMEHTS  HERE  INVALID 

2  IF  L-MOMEHTS  HERE  HOT  COHSISTEHT  WITH  THE 

RESTRICTIOHS  -l.LT.BETA.LT.l,  -1 .LT. DELTA. LT. 1 

3  IF  ESTIMATES  COULD  HOT  BE  FOUMD  ( ITERATION  FAILED 

TO  CONVERGE) 


***************************************************************************** 


SUBROUTINE  PELTU4(XM0M,  PARA,  GLD,  IFAIL) 

******  DECLARE  VARIABLES  AID  INITIALIZE  PARAMETERS 
IMPLICIT  REAL*8  (A-H.O-Z) 

REAL*8  XMOM(4) , PARA (5) ,  GLD (4) 

DATA  EPS/ID-6/ .MAXIT/20/ 

*  EPS  IS  TOLERANCE  OF  APPROXIMATION 

*  MAXIT  IS  MAX.  HUMBER  OF  ITERATIONS 

*  LOAD  LOCAL  VARIABLES  WITH  INPUT 

T3=XMQM(3) 

T4=XM0M(4) 

*  INITIALIZE  OUTPUT  ARRAY  WITH  ZEROES 

DO  10  1=1,5 
10  PARA(I)=ODO 

******  test  FOR  FEASIBILITY  AND  ERRORS 

*  CHECK  FOR  VARIANCE  LESS  THAN  OR  EQUAL  TO  ZERO 

IF(XM0M(2) .LE.ODO)GOTO  1000 

IF (DABS (T3) . GE . IDO . OR . DABS (T4) . GE . IDO) GOTO  1000 
IF (T4 . LT . 1 . 25*T3*T3-0 . 25) GOTO  1000 
IF (T4 . LT . ODO) GOTO  1010 
X=DABS(T3) 

IF (X* (5*X-T4) * (3*X*X-X*T4+2*T4*T4) .GT. (5*T4-X)*(X+T4)**2)G0T0  1010 

*  SET  STARTING  VALUES  FOR  N-R  ITERATION 

*  (ALGORITHM  CONVERGES  FAST,  SO  EVEN  A  CRUDE  GUESS  IS  OK) 

*  INITIAL  APPROXIMATIONS  OF  BETA  AND  DELTA 


B=ODO 


D=ODO 

******  NEWTON-RAPHSON  ITERATION 
DO  100  IT-1,  MAXIT 


N-2 


******  CALCULATE  CURREHT  ESTIMATES  OF  TAU-3  AID  TAU-4 
B1-1D0/(1D0+B) 

B2-1D0/C2D0+B) 

B3-1D0/C3D0+B) 

B4-1D0/(4D0+B) 

D1-1D0/(1D0+D) 

D2-1D0/(2D0+D) 

D3*1D0/(3D0+D) 

D4-1D0/(4D0+D) 

X2-B1*B2 

X3«*X2* ( 1DO-B) *B3 

X4»X3*(2D0-B)*B4 

Y2»D1*D2 

Y3»Y2*(1D0-D)*D3 

Y4=»Y3*  (2DO-D)  *D4 

AL2*  X2+Y2 

AL3— X3+Y3 

AL4=  X4+Y4 

TAU3-AL3/AL2 

TAU4=AL4/AL2 

******  JEST  FOR  COHVERGEMCE 

*  DIFFEREICE  BETWEEN  INPUT  AID  ESTIMATED  TAU-3  ft  TAU-4  VALUES 

E1=TAU3-T3 

E2=TAU4-T4 

IF (DABS (El) . LT. EPS. AID. DABS (E2) .LT.EPS)GOTO  120 

*  -  NOT  CONVERGED:  CALCULATE  NEXT  STEP 

*  -  NOTATION: 

*  DL2B  -  DERIVATIVE  OF  LAMBDA-2  W.R.T.  B 

*  DT3B  -  DERIVATIVE  OF  TAU-3  W.R.T.  B 

*  G..  -  MATRIX  OF  DERIVATIVES 

*  H. .  -  INVERSE  OF  G 

*  DEL.  -  STEPLENGTH 

BB1-1D0/ (1D0-B) 

BB2=1D0/ (2D0-B) 

DD1»1D0/(1D0-D) 

DD2=1D0/(2D0-D) 

DL2B=-X2*(B1+B2) 

DL3B=  X3*(B1+B2+B3+BB1) 

DL4B=-X4*(B1+B2+B3+B4+BB1+BB2) 

DL2D=-Y2* (D1+D2) 

DL3D=-Y3* (D1+D2+D3+DD1 ) 

DL4D=-Y4* (D1+D2+D3+D4+DD1+DD2) 
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DT3B= (DL3B-TAU3*DL2B) /AL2 
DT3D= (DL3D-TAU3*DL2D) /AL2 
DT4B= (DL4B-TAU4*DL2B) / AL2 
DT4D» (DL4D-TAU4*DL2D) /AL2 

*  MATRIX  6 

G11=DT3B 

G12-DT3D 

G21=»DT4B 

G22-DT4D 

*  DETERMINANT  OF  G 

DET=G11*G22-G12*G21 

*  MATRIX  H  (  INVERSE  OF  MATRIX  G  ) 

Hll=  G22/DET 
H12=-G12/DET 
H21=-G21/DET 
H22-  Gll/DET 

*  PRODUCT  OF  VECTOR  E  AND  MATRIX  H 

DEL1=E1*H11+E2*H12 

DEL2=E1*H21+E2*H22 

*  -  TAKE  THE  STEP 
B=B-DEL1 
D=D-DEL2 

*  -  IF  GONE  TOO  FAR,  REDUCE  STEPLENGTH 

*  WITH  UP  TO  TEN  HALF-STEPS  BACKWARDS 

IF (DABS (B) . LE . IDO . AND . DABS (D) . LE . IDO) GOTO  100 
DO  50  1=  1,10 

DEL 1=0 . 5* DELI 
DEL2=0 . 5*DEL2 
B=B+DEL1 
D=D+DEL2 

IF (DABS (B) . LE . IDO . AND . DABS (D) . LE . IDO) GOTO  100 
50  CONTINUE 

*  END  OF  N-R  ITERATION 
100  CONTINUE 

*  NOT  CONVERGED 
IFAIL=3 

RETURN 

*  CONVERGED 
120  IFAIL=0 

PARA(5)=D 

PARA (4) =XMOM( 2) / ( (B*X2 ) + (D*Y2) ) 

PARA(3)=B 
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PARA(2)*PARA(4) 

PARA(1)»XM0M(1)+PARA(2)*(D1-B1) 

******  COMPUTE  BERGEVIN ’ S  GLO  PARAMETERS  FROM  HOSKING’S  PARAMETERS 
GLD(1)SPARA(1) 

GLD(2)=1 .0D0/PARA(2) 

GLD (3) =PARA (3) 

GLD(4)SSPARA(5) 

RETURN 

******  error  traps 

1000  IFAIL=1 
RETURN 
1010  IFAIL*2 
RETURN 
END 
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Appendix  0.  FORTRAN  Programs:  L-Moment  to  GLD  Parameter 

Computer 

********4^**********************:|^******************************************* 

*  FORTRAM  77  CODE  USED  IN  AFIT/ENS  MASTER'S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  POWELL 

*  SUBROUTINE  PURPOSE:  COMPUTE  FOUR  GLD  PARAMETERS  FROM  FIRST  FOUR  L -MOMENTS 

*  USING  POVELL'S  ALGORITHM 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  searches  for  the  four  parameters  of  the 

*  Generalized  Lambda  Distribution  corresponding  to 

*  specified  values  for  the  resulting  1-moments 

*  of  the  distribution. 

*  AUTHOR(S) : 

*  E.  F.  MYKYTKA 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  46433,  USA 

*  DATE:  MARCH  1993 

*  INPUT  PARAMETERS: 

*  LMOM  -  DOUBLE  PRECISION  ARRAY  of  four  L-moments 

*  OUTPUT  PARAMETERS: 

*  LAMBDA  -  DOUBLE  PRECISION  ARRAY  of  four  GLD  parameters 

SUBROUTINE  POWELL  (  LMOM,  LAMBDA  ) 

C  DECLARE  VARIABLES 

IMPLICIT  REAL *8  (A-H,  0-Z) 

REAL+8  LAMBDA (4) ,  MEAN,  LM0M(4) 

DIMENSION  X(2),  E(2) ,  W(10) ,  START (2) 

COMMON  ALPHA3 ,  A3,  ALPHA4 ,  A4,  MEAN,  VAR 
C  E  —  Vector  ofvalues  used  to  determine 

C  convergence  of  Powell’s  algorithm 

C  ESCALE  —  Step-size  multiplier 

C  MAXIT  —  Maximum  number  of  iterations  for 

C  Powell ’ s  algorithm 

C  START  —  Vector  of  starting  values  for  the  lambda 

C  parameters  used  by  search  routine 

C  INITIALIZE  TOLERANCE  VALUE 
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TOL  -  2.D-4 

C  MAXIMUM  HUMBER  OF  ITERATIOHS 
MAXIT  -  500 

C  STEP-SIZE  MULTIPLIER 
ESCALE  »  1000 

C  COHVERGEHCE  RESOLUTIOH  VALUES 
E(l)  ■  0.000010D0 
E(2)  -  O.OOOOIODO 

C  LOAD  LOCAL  VARIABLES  WITH  INPUT  PARAMETERS 
MEAN  *  LMOM(l) 

VAR  =  LM0M(2) 

A3  =  LM0M(3) 

A4  *  LM0M(4) 

C  STARTING  POINTS  FOR  SEARCH  ROUTINE 
START(l)  =  0.050D0 
START(2)  =  0 . 05000 

C  INITIALIZATION  FOR  POWELL ’ S  ALGORITHM 
X(l)  -  START(l) 

X(2)  »  START (2) 

NI  *  5 
HO  -  6 

C  CALL  SUBROUTINE  BOTM  TO  PERFORM  POWELL  *  S  ALGORITHM 
CALL  BOTM(X,E,2,EF,ESCALE,3,MAXIT,W,NI,NO,10) 

C  INSURE  THAT  ALPHA3  AND  ALPHA4  CORRESPOND 

C  TO  OPTIMAL  VALUES  FOR  LAMBDA (3)  AND  LAMBDA (4) 

C  RETURNED  FROM  SUBROUTINE  BOTM 
CALL  CALCFX(2,  X,  EF) 

LAMBDA (3)  =  X(l) 

LAMBDA (4)  »  X(2) 

C  CALCULATE  LAMBDA(l)  AND  LAMBDA (2)  FROM  LAMBDA (3) 

C  AND  LAMBDA (4) 

CALL  FN2 (LAMBDA) 

RETURN 

END 

♦♦♦*****♦**♦#*♦♦#♦**#* #♦♦♦**♦♦♦*** *♦***♦#*****♦♦*****♦**♦*#♦****#♦♦♦♦♦*#♦**** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  CALCF7 

*  SUBROUTINE  PURPOSE:  CALCULATE  THE  OBJECTIVE  FUNCTION  FOR  POWELL’S 

*  ALGORITHM 
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*  SUBROUTIVE  DESCRIPTION : 

*  Evaluates  the  objective  function  to  be  minimized, 

*  which  is  the  sum  of  the  squared  differences  between 

*  the  desired  and  calculated  values  of  skewness  and 

*  kurtosis 

*  AUTHOR(S) : 

*  E.  F.  MYKYTKA 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  HRIGHT-PATTERSOH  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  MARCH  1993 

*  INPUT  PARAMETERS: 

*  L  -  DOUBLE  PRECISION  ARRAY  two  variables  of  function 

*  OUTPUT  PARAMETERS: 

*  FN  -  DOUBLE  PRECISION  function  value 

*  N  -  INTEGER  count  of  subroutine  calls 
******************************* *** ******************************************* 

SUBROUTINE  CALCFX  (  N,  L,  FN  ) 

C  DECLARE  VARIABLES 

IMPLICIT  REAL *8  (A-H,  0-Z) 

REAL*8  L(2) ,  LAMBDA (4) ,  MEAN 
COMMON  ALPHA3,  A3,  ALPHA4,  A4,  MEAN,  VAR 
C  INSURE  THAT  LAMBDA (3)  AND  LAMBDA (4)  HAVE  SAME  SIGN 
LAMBDA (3)  =  L(l) 

LAMBDA (4)  -  L(2) 

IF  (LAMBDA(3)*LAMBDA(4) .LT.ODO)  GOTO  999 
C  COMPUTE  OBJECTIVE  FUNCTION 
B=L(1) 

D=L(2) 

B1=1D0/ (1DO+B) 

B2=1D0/ (2DO+B) 

B3=1D0/ (3D0+B) 

B4=1D0/ (4D0+B) 

D1=1D0/ (1DO+D) 

D2=1D0/ (2D0+D) 

D3=1D0/ (3D0+D) 

D4=1D0/ (4D0+D) 

X2=B*B1*B2 
X3=X2* ( 1DO-B) *B3 
X4SSX3*  (2DO-B)  *B4 
Y2=D*D1*D2 
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Y3-Y2*(1D0-D)*D3 
Y4»Y3* (2D0-D) *D4 
AL2-  X2+Y2 
AL3—X3+Y3 
AL4-  X4+Y4 
ALPHA3-AL3/AL2 
ALPHA4-AL4/AL2 

C  DIFFERENCE  BETWEEN  DESIRED  AND  CALCULATED  VALUES 
DIFF3  -  ALPHAS  -  A3 
DIFF4  -  ALPHA4  -  A4 

C  FUNCTION  VALUE  IS  SQUARED  DIFFERENCE  BETWEEN 

C  CALCULATED  AND  DESIRED  ALPHA3  AND  ALPHA4 
FN  *  DIFF3*DIFF3  +  DIFF4*DIFF4 
RETURN 

999  FH  =  10. DO 
RETURN 
END 

#****#♦♦♦##♦#♦**#*+**+****++****#♦♦**♦*************+****+**+***#+******#♦**** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE  :FN2 

*  SUBROUTINE  PURPOSE:  CALCULATE  LAMBDA-1  AND  LAMBDA-2  AS  FUNCTIONS  OF 

*  LAMBDA-3  AND  LAMBDA-4 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  receives  an  incomplete  array  of  Lambda  Parameters  and 

*  completes  the  array  as  a  function  of  the  parameters  provided. 

*  AUTHOR(S) : 

*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT -PATTERS ON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  E.  F.  MYKYTKA 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  OCTOBER  1993 

*  INPUT  PARAMETERS : 

*  LAMBDA  -  DOUBLE  PRECISION  ARRAY  of  four  Lambda  Parameters 
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*  OUTPUT  PARAMETERS : 

*  LAMBDA  -  DOUBLE  PRECIS IOH  ARRAY  of  foui  Lambda  Parameters 
***************************************************************************** 

SUBROUTINE  FV2  (  LAMBDA  ) 

C  DECLARE  VARIABLES 

IMPLICIT  REAL*8  (A-H,  0-Z) 

REAL*8  MEAN,  LAMBDA (4) 

COMMON  ALPHA3,  A3,  ALPHA4,  A4,  MEAN,  VAR 
C  COMPUTE  LAMBDA-1  AND  LAMBDA-2 
B  =  LAMBDA (3) 

D  *  LAMBDA (4) 

B1»1D0/(1D0+B) 

B2=1D0/(2D0+B) 

D1*1D0/(1D0+D) 

D2»1D0/(2D0+D) 

X2*B1*B2 

Y2=SD1*D2 

C  LOAD  OUTPUT  PARAMETERS  AND  RETURN 
LAMBDA (2) =VAR/ ( (B*X2) + (D*Y2) ) 

LAMBDA (1 ) =MEAN+LAMBDA (2) * (Dl-Bl) 

LAMBDA (2) =1 . ODC/LAMBDA (2) 

RETURN 

END 
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Appendix  P.  FORTRAN  Program:  PowelVs  Algorithm  for  function 

minimization 

*ii^**********************im************************************************** 

*  FORTRAH  77  CODE  USED  IN  AFIT/EHS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMEHTS  TO  FIT  THE  GEIERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTIME  TITLE:  BOTM 

*  SUBROUTIME  PURPOSE:  MIMIMIZE  OBJECTIVE  FUMCTIOH  USIMG  POWELL’S  ALGORITHM 

*  SUBROUTIME  DESCRIPTION: 

*  Subroutine  performs  Powell’s  Algorithm  for  function 

*  minimization,  developed  by  M.  J.  D.  Powell 

*  Source:  Keuster,  J.  L.  and  Mize,  J.  H., 

*  OPTIMIZATION  TECHNIQUES  WITH  FORTRAN, 

*  New  York:  McGraw-Hill,  1973 

*  AUTHOR(S) : 

*  E.  F.  MYKYTKA 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  HRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  MARCH  1993 

*  INPUT / OUTPUT  PARAMETERS: 

*  X  -  DOUBLE  PRECISION  ARRAY  OF  c  ,  function  variables 

*  E  -  DOUBLE  PRECISION  ARRAY  OF  S  2,  convergence  variables 

*  N  -  INTEGER  VARIABLE,  2 

*  EF  -  DOUBLE  PRECISION  VARIABLE 

*  ESCALE  -  DOUBLE  PRECISION  VARIABLE,  step-size  multiplier 

*  IPRIHT  -  INTEGER  VARIABLE 

*  MAXIT  -  INTEGER  VARIABLE,  maximum  number  of  iterations 

*  W  -  DOUBLE  PRECISION  ARRAY  OF  SIZE  NH 

*  NI  -  INTEGER  VARIABLE,  5 

*  NO  -  INTEGER  VARIABLE,  6 

*  NH  -  INTEGER  VARIABLE,  10 

**************1^************************************************************* 
SUBROUTINE  BOTM  (X,  E,  N,  EF,  ESCALE,  IPRINT, 

+MAXIT,  H,  NI,  NO,  NH) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  X(N) ,  H(NH) ,  E(N) 

DDMAG=0.1D0*ESCALE 
SCER=0 . 05D0 /ESCALE 
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JJ-i*(S+l) 

JJJ-JJ+I 

K-M+l 

MFCC-1 

IID-1 

m-i 

DO  4  1*1,1 
H(I)*ESCALE 
DO  4  J*1,I 
H(K)*O.DO 
IF(I-J)4,3,4 

3  H(K)*DABS  (E(I)) 

4  K=K+1 
ITEROl 
IS GRAD* 2 

CALL  CALCFX(N,X,F) 

FKEEP=2 . DO*DABS  (F) 

5  ITOHE=l 
FP*F 
SUM-O. DO 
IXP=JJ 

DO  6  1*1 ,M 
IXP-IXP+1 

6  W(IXP)*X(I) 

IDIRH=M+1 

ILINE-l 

7  DMAX* V ( ILIHE) 

DACC=DMAX* S CER 

DMAG-DMIN1  (DDMAG , 0 . 1D0*DMAX) 
DMAG-DMAX1 (DMAG , 20 . DO*DACC) 
DDMAX-10 . DO*DMAG 
GO  TO  (70,70,71) ,ITONE 
70  DL-O.DO 
D-DMAG 
FPREV-F 
IS=5 

FA*FPREV 

DA-DL 

8  DD=D-DL 
DL=D 

58  K=IDIRN 
DO  9  1*1, H 
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X(I)«X(I)+DD*H(K) 

9  K-K+l 

CALL  CALCFX(M,X,F) 

HFCC-RFCC+1 

GO  TO  (10, 11, 12, 13, 14, 96), IS 

14  IF (F-FA) 15,16, 24 

16  IF  (DABS  (D)-DMAX)  17,17,18 

17  D-D+D 
GO  TO  8 

18  WRITE(*,019) 

19  FORMAT  (5X , 38HMAXIMUM  CHARGE  DOES  ROT  ALTER  FURCTIOR) 
GO  TO  20 

15  FB«F 
DB=D 

GO  TO  21 

24  FB=FA 
DB-DA 
FA-F 
DA-D 

21  GO  TO  (83,23) ,ISGRAD 
23  D-DB+DB-DA 
IS=1 
GO  TO  8 

83  D=0 . 6D0* (DA+DB- (FA-FB) / (DA-DB) ) 

IS =4 

IF ( (DA-D) * (D-DB) ) 25 , 8 , 8 

25  IS-1 

IF (DABS  (D-DB) -DDMAX) 8,8,26 

26  D-DB+DSIGR  (DDMAX, DB-DA) 

IS*1 

DDMAX=DDMAI+DDMAX 

DDMAG=DDMAG+DDMAG 

IF  (DDMAG.GE.l.D+60)  DDMAG  =  l.D+60 
IF (DDMAX-DMAX) 8,8,27 

27  DDMAX=DMAX 
GO  TO  8 

13  IF (F-FA) 28, 23, 23 

28  FC=FB 
DC=DB 

29  FB=F 
DB=D 

GO  TO  30 
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12  IF (F-FB) 28,28, 31 

31  FA-F 
DA-D 

60  TO  30 

11  IF (F-FB) 32, 10, 10 

32  FA-FB 
DA-DB 

GO  TO  29 
71  DL-l.DO 
DDMAX-5 . 

FA-FP 
DA— l.DO 
FB-FHOLD 
DB-O.DO 
D=1 . DO 
10  FOF 
DC-D 

30  A- (DB-DC) * (FA-FC) 

B- (DC-DA) * (FB-FC) 

IF ( (A+B) * (DA-DC) ) 33 , 33 , 34 

33  FA-FB 
DA-DB 
FB-FC 
DB-DC 

GO  TO  26 

34  D-0 . 6D0* (A* (DB+DC) +B* (DA+DC) ) / (A+B) 

DI-DB 

FI-FB 

IF(FB-FC)44,44,43 

43  DI-DC 
FI-FC 

44  GO  TO  (86,86,85), ITOHE 

85  ITOME-2 
GO  TO  45 

86  IF  (DABS  (D-DI)-DACC)  41,41,93 

93  IF  (DABS  (D-DI) -0 . 03D0*DABS  (D))  41,41,45 

45  IF  ( (DA-DC) *(DC-D))  47,46,46 

46  FA-FB 
DA-DB 
FB-FC 
DB-DC 

GO  TO  25 
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47  IS-2 

IF  ( (DB-D) * (D-DC) )  48,8,8 

48  IS-3 

60  TO  8 
41  F=FI 
D-DI-DL 

DD*DSQRT  ( (DC-DB) * (DC-DA) * (DA-DB) / (A+B) ) 

DO  49  1=1, H 

X(I)*X(I)+D*W(IDIRH) 

V(IDIRH) *DD*V ( ID IRH) 

49  IDIRH=IDIRH+1 

W ( ILIHE) *H ( ILIHE) /DD 
ILINE-ILIHE+1 
IF(IPRIHT-l) 61 ,60,51 

50  corriHUE 

GO  TO  (51 , 53 , 53) , IPRIMT 

51  60  TO  (55,38) ,ITOHE 

65  IF  (FPREV-F-SUM)  94,95,95 

95  SUM*FPREV-F 
JIL*ILIHE 

94  IF  (IDIRM-JJ)  7,7,84 
84  60  TO  (92,72) ,IHD 
92  FHOLD=F 
IS=6 
IXP*JJ 
DO  59  1=1, N 
IXP=IXP+1 

59  H(IXP)=X(I)-W(IXP) 

DD=1 . DO 

60  TO  58 

96  GO  TO  (112,87) ,IHD 
112  IF(FP-F)  37,37,91 

91  D=2 . DO* (FP+F-2 .DO*FHOLD) / (FP-F) **2 
IF  (D* (FP-FHOLD-SUM) **2-SUM)  87,37,37 
87  J=JIL*H+1 

IF  (J-JJ)  60,60,61 

60  DO  62  I=J,JJ 
K=I-N 

62  W(K)=W(I) 

DO  97  I=JIL,N 

97  W(I-1)=W(I) 

61  IDIRN=IDIRN-N 
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ITONE-3 

X-IDIRH 

IXP-JJ 

AAA-0 . DO 

DO  67  1-1,1 

IXP-IIP+1 

W(K)-W(IXP) 

IF  (AAA-DABS  (W(K)/E(I)))  66,67,67 

66  AAA-DABS  (W(K)/E(I)) 

67  K-K+l 
DDHAG-l.DO 

W  (N) -ESC ALE/ AAA 
ILINE-N 
GO  TO  7 

37  IXP-JJ 
AAA-0 . DO 
F-FHOLD 

DO  99  1-1 ,H 

IXP-IXP+1 

X(I)-X(I)-W(IXP) 

IF(AAA*DABS  (E(I))-DABS  (W(IXP)))  98,99,99 

98  AAA-DABS  (W(IXP)/E(I)) 

99  CONTINUE 
GO  TO  72 

38  AAA-AAA* (1 . DO+DI) 

GO  TO  (72,106) ,I1D 

72  IF ( IPRIMT-2) 53,50, 50 
53  GO  TO  (109,88) ,IHD 
109  IF  (AAA  -  0 . IDO)  20,20,76 
76  IF (F-FP) 35,78, 78 
78  WRITE (*,80) 

80  F0RMAT(5X,31HACCURACY  LIMITED  BY  ERRORS  IN  F) 
GO  TO  20 

88  IHD-1 

35  DDMAG-0 . 4D0*DSQRT (DABS (FP-F) ) 

IF  (DDMAG . GE . 1 . D+60)  DDMAG  =  l.D+60 
ISGRAD-1 

108  ITERC-ITERC+1 

IF ( ITERC-MAXIT) 5,5,81 

81  WRITE (*,82)  MAXIT 

82  FORMAT ( 15, 29H  ITERATIONS  COMPLETED  BY  BOTM) 

IF (F-FKEEP) 20,20, 110 
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110  F-FKEEP 

DO  111  1-1,1 
JJJ-JJJ+1 

111  X(I)-W(JJJ) 

GO  TO  20 

101  JIL-1 
FP-FKEEP 

IF (F-FKEEP ) 105 , 78 , 104 

104  JIL-2 
FP-F 
F-FKEEP 

105  IXP-JJ 

DO  113  1-1,1 

IXP-IXP+1 

K-IXP+M 

GO  TO  (114,115) , JIL 
114  V(IXP)-W(K) 

GO  TO  113 
116  W(IXP)=X(I) 

X(I)-H(K) 

113  CONTINUE 
JIL-2 
GO  TO  92 

106  IF (AAA-0 . IDO)  20,20,107 
20  EF-F 

RETURN 

107  INM-1 

GO  TO  36 
END 
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Appendix  Q.  FORTRAN  Programs:  Report  File  Generator 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IM  AFIT/ENS  MASTER'S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GEHERALIZED  LAMBDA  DISTRIBUTIOI 

*  TO  SAMPLE  DATA" 

*  SUBROUTIHE  TITLE  :  REPORT 

*  SUBROUTIHE  PURPOSE:  GENERATE  REPORT  OUTPUT  FILE 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  collects  all  the  information  generated  by  the  program 

*  and  outputs  it  to  a  text  fil~  in  a  readable  and  labeled  format.  The 

*  subroutine  prompt s  the  user  for  a  filename.  Output  includes: 

*  number  of  elements  in  sample 

*  First  four  L -moments  of  sample 

*  Hosking' 8  Lambda  Parameters  from  Heston-Raphson  Search 

*  Ramberg’s  Lambda  Parameters  derived  from  Hosking’s  Parameters 

*  Ramberg’s  Lambda  Parameters  from  Posell’s  Algorithm 

*  Re-computed  L-moments  from  H-R  search 

*  Re-computed  L-moments  from  Posell’s  algorithm 

*  AUTHOR(S) : 

*  R.  B.  MOHAN,  CAPTA TH ,  USAF 

*  DEPARTMENT  OF  OPER  '.riONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSOH  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  OCTOBER  1993 

*  INPUT  PARAMETERS: 

*  HUM:  INTEGER  PARAMETER  number  of  points  in  set 

*  I,  II:  INTEGER  PARAMETER  failure  codes 

*  MOMENT:  DOUBLE  PRECISION  PARAMETER  data  array  -  L-moments 

*  LAMBDA:  DOUBLE  PRECISION  PARAMETER  data  array  -  Hosking’s  parameters 

*  NEVA  DOUBLE  PRECISION  PARAMETER  data  array  -  N-R  L-moments 

*  NEVB:  DOUBLE  PRECISION  PARAMETER  data  array  -  Powell  L-moments 

*  GLD:  DOUBLE  PRECISION  PARAMETER  data  array  -  Ramberg’s  parameters 

*  POWGLD:  DOUBLE  PRECISION  PARAMETER  data  array  -  Ramberg’s  parameters 

*  OUTPUT  PARAMETERS :  NONE. 

* 

***************************************************************************** 
SUBROUTINE  REPORT  (NUM, MOMENT .LAMBDA, GLD, POWGLD, NEVA, NEVB, I) 

******  DECLARE  LOCAL  VARIABLES 

DOUBLE  PRECISION  MOMENT  (4) ,  LAMBDA  (5) ,  NEVA  (4) ,  GLD  (4) 

DOUBLE  PRECISION  POWGLD  (4) ,  NEVB  (4) 
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INTEGER  NUM,  I,  II 
CHARACTER*  12  NAME 

******  PROMPT  USER  FOR  NAME  OF  FILE 

PRINT  *,  ’ENTER  THE  NAME  OF  THE  OUTPUT  FILE:’ 

READ  (*,  ’ (A12) ’  )  NAME 

******  open  OUTPUT  FILE  AND  DOWNLOAD  FORMATTED  DATA 

OPEN  (UNIT=30 ,  FILENAME,  STATUS=’NEW’ ,  IQSTAT=II,  ERR=200) 


******  format 


10 

20 

30 


FORMAT  (10X) 


FORMAT 

FORMAT 


( 


170 

180 

210 

220 

230 

240 

250 

260 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 


NUMBER  OF  ELEMENTS  IN  SAMPLE  DATA:  ’ 
FIRST  FOUR  MOMENTS  OF  SAMPLE  DATA’) 


14) 


40 

FORMAT 

( 

’L-l :  ’ 

,  G13.6) 

50 

FORMAT 

( 

’L-2:  ’ 

,  G13.6) 

60 

FORMAT 

( 

’T-3:  ’ 

,  G13.6) 

70 

FORMAT 

( 

’T-4:  ’ 

,  G13.6) 

80 

FORMAT 

( 

’LAMBDA  DISTRIBUTION  PARAMETERS 

90 

FORMAT 

( 

’XI: 

’,  G13.6) 

100 

FORMAT 

( 

’ALPHA: 

’,  G13.6) 

110 

FORMAT 

( 

’BETA: 

’,  G13.6) 

120 

FORMAT 

( 

’GAMMA: 

’,  G13.6) 

130 

FORMAT 

( 

’DELTA: 

’,  G13.6) 

140 

FORMAT 

( 

’L-MOMENTS  ARE  INVALID’) 

150 

FORMAT 

( 

’L-MOMENTS  ARE  INCONSISTENT’) 

160 

FORMAT 

( 

’FAILURE  TO  CONVERGE’) 

RE-COMPUTED  L-MOMENTS  (Newton-Raphuon) ’ ) 
GLD  DISTRIBUTION  PARAMETERS  (Bergevin) ’) 
Lambda-1:  ’,  G13.6) 

Lambda-2:  ’,  G13.6) 

Lambda-3:  ’,  G13.6) 

Lambda-4:  ’,  G13.6) 

GLD  DISTRIBUTION  PARAMETERS  (Powell)’) 
RE-COMPUTED  L-MOMENTS  (Powell) ’ ) 


******  write  to  output  file 
WRITE  (30,10) 

WRITE  (30,10) 

WRITE  (30,20)  NUM 
WRITE  (30,10) 

WRITE  (30,30) 

WRITE  (30,10) 

WRITE  (30,40)  MOMENT(l) 
WRITE  (30,10) 

WRITE  (30,50)  MOMENT (2) 
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WRITE  (30,10) 

WRITE  (30,60)  MOMENT (3) 
WRITE  (30,10) 

WRITE  (30,70)  MOMENT (4) 
WRITE  (30,10) 

IF  (I  .GT.  0)  GOTO  190 
WRITE  (30,80) 

WRITE  (30,10) 

WRITE  (30,90)  LAMBDA (1) 
WRITE  (30,10) 

WRITE  (30,100)  LAMBDA (2) 
WRITE  (30,10) 

WRITE  (30,110)  LAMBDA (3) 
WRITE  (30,10) 

WRITE  (30,120)  LAMBDA (4) 
WRITE  (30,10) 

WRITE  (30,130)  LAMBDA(5) 
WRITE  (30,10) 

WRITE  (30,180) 

WRITE  (30,10) 

WRITE  (30,210)  GLD(l) 
WRITE  (30,10) 

WRITE  (30,220)  GLD(2) 
WRITE  (30,10) 

WRITE  (30,230)  GLD(3) 
WRITE  (30,10) 

WRITE  (30,240)  GLD(4) 
WRITE  (30,10) 

WRITE  (30,250) 

WRITE  (30,10) 

WRITE  (30,210)  POWGLD(l) 
WRITE  (30,10) 

WRITE  (30,220)  P0WGLD(2) 
WRITE  (30,10) 

WRITE  (30,230)  P0WGLD(3) 
WRITE  (30,10) 

WRITE  (30,240)  P0WGLD(4) 
WRITE  (30,10) 

WRITE  (30,170) 

WRITE  (30,10) 

WRITE  (30,40)  NEWA(l) 
WRITE  (30,10) 
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WRITE  (30,50)  NEVA (2) 

WRITE  (30,10) 

WRITE  (30,60)  MEWA(3) 

WRITE  (30,10) 

WRITE  (30,70)  NEWA (4) 

WRITE  (30,10) 

WRITE  (30,260) 

WRITE  (30,10) 

WRITE  (30,40)  NEWB(l) 

WRITE  (30,10) 

WRITE  (30,50)  HEWB(2) 

WRITE  (30,10) 

WRITE  (30,60)  NEWB(3) 

WRITE  (30,10) 

WRITE  (30,70)  NEWB (4) 

WRITE  (30,10) 

WRITE  (30,10) 

190  IF  (I.Eq.l)  WRITE  (30,140) 

IF  (I.Eq.2)  WRITE  (30,150) 

IF  (I.Eq.3)  WRITE  (30,160) 

******  close  output  file 
CLOSE  (30) 

******  RETURN  TO  MAIN  PROGRAM 
RETURN 

******  ERROR  TRAP 

200  PRINT  *,  1  CANNOT  OPEN  FILE  \  NAME,  »,  ERROR=  *,  II 
RETURN 
END 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER'S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  RECOMP 

*  SUBROUTINE  PURPOSE:  RE-COMPUTES  L-MOMENTS  FROM  LAMBDA  PARAMETERS  TO 

*  DOUBLE-CHECK  ACCURACY  OF  SEARCH  ROUTINES 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  receives  the  four  lambda  parameters  from  a  search  routine 

*  and  runs  them  back  through  the  defining  equations  for  the  L-moments  to 

*  recompute  the  L-moments.  This  acts  as  a  check  of  the  accuracy  of  the 

*  search  routines . 

*  AUTHOR(S) : 
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*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  HRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  OCTOBER  1993 

*  INPUT  PARAMETERS : 

*  LAMBDA  -  Double  Precision  Array  of  four  lambda  parameters 

*  OUTPUT  PARAMETERS: 

*  NEVMOM  -  Double  Precision  Array  of  four  L-moments 

***************************************************************************** 

SUBROUTINE  RECOMP  (  LAMBDA,  NEVMOM  ) 

******  DECLARE  LOCAL  VARIABLES 

DOUBLE  PRECISION  LAMBDA  (  1:5  ),  NEVMOM  (  1:4  ) 

DOUBLE  PRECISION  LI,  L2,  L3,  L4,  T3,  T4,  X,  A,  B,  C,  D,  PI,  P2 

******  INITIALIZE  LOCAL  VARIABLES 
X  =  LAMBDA  (1) 

A  »  LAMBDA  (2) 

B  =  LAMBDA  (3) 

C  =  LAMBDA  (4) 

D  =  LAMBDA  (5) 

******  CALCULATE  L-MOMENTS  FROM  LAMBDA  PARAMETERS 
PI  =  (A  /  (1+B) ) 

P2  =  (C  /  (1+D)) 

LI  ■  X  +  PI  -  P2 
PI  =  PI  *  (B  /  (2+B)) 

P2  -  P2  *  (D  /  (2+D)) 

L2  =  PI  +  P2 

PI  -  PI  *  ( (B— 1 )  /  (3+B) ) 

P2  =  P2  *  ((D-l)  /  (3+D)) 

T3  -  (  PI  -  P2  )  /  L2 
PI  =  PI  *  C(B-2)  /  (4+B)) 

P2  *  P2  *  ( (D-2)  /  (4+D)) 

T4  =  (  PI  +  P2  )  /  L2 

******  LOAD  OUTPUT  PARAMETERS 
NEVMOM  (1)  =  LI 
NEVMOM  (2)  -  L2 
NEVMOM  (3)  =  T3 
NEVMOM  (4)  =  T4 

******  RETURN  TO  MAIN  PROGRAM 
RETURN 
END 


Q-5 


Appendix  R.  FORTRAN  Program:  PDF  Plot  File  Generator 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IM  AFIT/ENS  MASTER'S  THESIS; 

*  "THE  USE  OF  L -MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  PLOTFILE 

*  SUBROUTINE  PURPOSE:  GENERATE  DATA  OUTPUT  FILE  TO  FEED  HATLAB  PLOTTING 

*  ROUTINE  TO  GENERATE  GRAPH 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  receives  the  four  GLD  parameters  and  generates  a  two- 

*  dimensional  array  of  plotting  pairs.  There  are  999  pairs  of  p  (0  to  1) 

*  vs.  R(p) ,  determined  by  the  GLD  percentile  function  and  the  input  GLD 

*  parameters.  The  subroutine  then  prompts  the  user  for  a  file  name  and 

*  outputs  the  array  to  the  file. 

*  AUTHOR(S) : 

*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  OCTOBER  1993 

*  INPUT  PARAMETERS: 

*  GLD  -  DOUBLE  PRECISION  ARRAY  of  four  Lambda  Parameters 

*  OUTPUT  PARAMETERS:  NONE. 

***************************************************************************** 
SUBROUTINE  PLOTFILE  (  GLD  ) 

******  DECLARE  LOCAL  VARIABLES 

DOUBLE  PRECISION  GLD  (  1:4  ) 

DOUBLE  PRECISION  PLOT  (1:999,  1:2  ) 

DOUBLE  PRECISION  LI,  L2,  L3,  L4,  P,  X,  FX 
INTEGER  NUM,  I,  II 
CHARACTER  ANSWER 
CHARACTER*  12  NAME 

******  INITIALIZE  LOCAL  VARIABLES 
P  =  ODO 
NUM  *  999 
LI  -  GLD  (1) 

L2  =  GLD  (2) 

L3  =  GLD  (3) 

L4  =  GLD  (4) 

******  COMPUTE  ARRAY  OF  PLOTTING  POINT  PAIRS 
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DO  300  I  *  1,  HUM 
P  =  P  +  0.001 ODO 

X  »  LI  +  ((P**L3)  -  ((l-P)**L4))/L2 

FX  -  L2/ ( (L3* (P** (L3-1) ) ) +(L4* ( (1-P) ** (L4-1) ) ) ) 

PLOT  (  I,  1  )  =  X 
PLOT  (  I,  2  )  =  FX 
300  COHTINUE 

******  PROMPT  USER  FOR  HAME  OF  FILE 

PRIHT  *,  'DO  YOU  WANT  TO  GENERATE  A  DATA  PLOTTING  FILE?' 
READ  *,  ANSWER 

IF  ( (ANSWER. NE. 'Y') .AND. (ANSWER. NE. 'y' ))  THEN 
GOTO  400 
END  IF 

PRINT  *,  'ENTER  THE  NAME  OF  THE  PLOT  FILE:' 

READ  (*,  ' (A12) '  )  NAME 

******  open  OUTPUT  FILE  AND  DOWNLOAD  FORMATTED  DATA 

OPEN  (UNIT=30 ,  FILE=NAME,  STATUS='NEW' ,  IOSTAT=II,  ERR=200) 
10  FORMAT  (F11.5,  2X,  F11.5) 

******  download  array 
DO  100  I  =  1,  NUM 

WRITE  (30,  10)  PLOT  (1,1),  PLOT  (1,2) 

100  CONTINUE 
******  close  output  file 
CLOSE  (30) 

******  RETURN  TO  MAIN  PROGRAM 
400  RETURN 
******  ERROR  TRAP 

200  PRINT  *,  'CANNOT  OPEN  FILE  ',  NAME,  ',  ERROR=  ',  II 
RETURN 
END 
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Appendix  S.  FORTRAN  Programs:  CDF  Plot  File  Generator 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-KOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE : CDFPLOT 

*  SUBROUTINE  PURPOSE:  BUILDS  PLOTFILE  FOR  CDF  OF  GLD  FROM  LAMBDA  PARAMETERS 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  builds  an  array  of  plotting  points  of  the  CDF  of  a 

*  theoretical  distribution  using  internal  lambda  parameters  and  an  array  of 

*  plotting  points  of  the  CDF  of  an  estimated  distribution  using  external 

*  lambda  parameters.  It  then  downloads  the  array  to  an  output  file.  It 

*  also  determines  Kolmogorov- Smirnov  Statistics  by  comparing  the  two 

*  GLD  CDF’s:  MIN,  MAX,  and  AVG  values  of  differences  between  theoretical  and 

*  empirical  CDF’s. 

*  AUTHOR(S) : 

*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

+  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSOH  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  NOVEMBER  1993 

*  INPUT  PARAMETERS: 

*  LAMBDA:  DOUBLE  PRECISION  PARAMETER  data  array 

*  OUTPUT  PARAMETERS:  NONE. 

***************************************************************************** 
SUBROUTINE  CDFPLOT  (  POWGLD  ) 

******  DECLARE  LOCAL  VARIABLES 

DOUBLE  PRECISION  POWGLD  (4) ,  THEOGLD  (4) 

DOUBLE  PRECISION  PLOT  (1:599,  1:3  ) 

DOUBLE  PRECISION  X,  FX1,  FX2,  DELTA,  MIN,  MAX,  AVG 

INTEGER  NUM,  I,  II 

CHARACTER  ANSWER 

CHARACTER* 12  NAME 

DOUBLE  PRECISION  CDF 

EXTERNAL  CDF 

******  INITIALIZE  LOCAL  VARIABLES  (normal  distribution) 

X  =  -3.0D0  *  first  plotting  pair 

HUM  =  599  *  number  of  plotting  pairs 

THEOGLD  (1)  =  O.ODO  *  theoretical  mean  (lambda  1) 

THEOGLD  (2)  =  0.19750D0  *  theoretical  lambda  2 
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THEOGLD  (3)  *  0.13490D0  *  theoretical  lambda  3 
THEOGLD  (4)  **  0.13490D0  *  theoretical  lambda  4 

HIM  -  1 .000 
MAI  *  0.0D0 
AVG  ■  0.000 

******  COMPUTE  ARRAY  OF  PLOTTING  POINT  PAIRS 
DO  300  1=1,  MUM 

X  «  X  +  0.010D0  *  increment  between  plotting  points 
FX1  *  CDF  (  THEOGLD,  X  )  *  theoretical  cdf  value 
FX2  =  CDF  (  POVGLD ,  X  )  *  empirical  cdf  value 
PLOT  (  I,  1  )  -  X 
PLOT  (I,  2  )  *  FX1 
PLOT  (  I,  3  )  =  FX2 

******  COMPUTE  KOLMOGOROV-SMIRHOV  STATISICS 

*  compute  the  difference  between  the  F(x)  values  for  theoretical  and 

*  empirical  distributions 

DELTA  =  DABS  (  FX1  -  FX2  ) 

*  search  for  the  min-imum  difference  between  the  distributions 

IF  (  DELTA  .LT.  MIN  )  MIN  =  DELTA 

*  search  for  the  lMT-iimim  difference  between  the  distributions 

IF  (  DELTA  .GT.  MAX  )  MAX  =  DELTA 

*  compute  the  average  of  all  the  differences  between  the  distributions 

AVG  *  AVG  +  DELTA 
300  CONTINUE 

AVG  *  AVG/HUM 

*  output  K-S  stats  to  screen 

PRINT  *,  »  ’ 

PRINT  *,  ’K-S  Statistics’ 

PRINT  *,  ’  ’ 

PRINT  *,  ’MINIMUM  DELTA-CDF:  ’,  MIN 
PRINT  *,  ’AVERAGE  DELTA-CDF:  ’,  AVG 
PRINT  *,  ’MAXIMUM  DELTA-CDF:  ’,  MAX 
PRINT  *,  ’  ’ 

******  PROMPT  USER  FOR  NAME  OF  FILE 

PRINT  *,  ’DO  YOU  WANT  TO  GENERATE  A  CDF  DATA  PLOTTING  FILE?’ 

READ  *,  ANSWER 

IF  ( (ANSWER. ME. ’Y’) .AND. (ANSWER. NE.’y’))  THEN 
GOTO  400 
END  IF 

PRINT  *,  ’ENTER  THE  NAME  OF  THE  CDF  PLOT  FILE:’ 

READ  (*,  ’ (A12) ’  )  NAME 

******  open  OUTPUT  FILE  AND  DOWNLOAD  FORMATTED  DATA 
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OPEN  (UNIT-30,  FILENAME,  STATUS-'NEW' ,  IOSTAT-II,  ERR=200) 
10  FORMAT  (F11.5,  2X,  F11.5,  2X,  FU.5) 

******  download  array 
DO  100  I  3  1,  MUM 

WRITE  (30,  10)  PLOT  (1,1),  PLOT  (1,2),  PLOT  (1,3) 

100  CONTINUE 
******  close  output  file 
CLOSE  (30) 

******  RETURN  TO  MAIN  PROGRAM 
400  RETURN 
******  ERROR  TRAP 

200  PRINT  *,  ’ CANNOT  OPEN  FILE  »,  NAME,  ',  ERROR=  \  II 
RETURN 
END 


***************************************************************************** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER'S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE: CDF 

*  SUBROUTINE  PURPOSE:  COMPUTE  F(X)  VALUES  OF  CDF  OF  GLD,  GIVEN  X  VALUES 

*  SUBROUTINE  DESCRIPTION: 

*  This  function  uses  a  secant  search  to  find  the  value  of  p  between  0  and 

*  1,  given  the  value  of  X.  It  then  finds  the  value  of  F(X)  using  the  derived 

*  value  of  p. 

*  AUTHOR(S) : 

*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  E.  F.  MYKYTKA 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  NOVEMBER  1993 

*  INPUT  PARAMETERS: 

*  LAMBDA:  DOUBLE  PRECISION  PARAMETER  data  array 

*  X:  DOUBLE  PRECISION  PARAMETER  independent  variable 

*  OUTPUT  PARAMETERS: 

*  CDF:  DOUBLE  PRECISION  PARAMETER  dependent  variable 
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***************************************************************************** 


FUMCTIOM  CDF  (  LAMBDA,  X  ) 

******  DECLARE  LOCAL  VARIABLES 

IMPLICIT  DOUBLE  PRECISIOM  (  A-H,  O-Z  ) 

DOUBLE  PRECISIOM  LAMBDA (4) 

DOUBLE  PRECISIOM  GLDQF 
EXTERMAL  GLDQF 

*  Check  to  see  if  the  range  is  bounded  or  not 

IF (LAMBDA (3). LT. O.O)  GOTO  20 
IF (LAMBDA (4). LT. 0.0)  GOTO  20 

*  Range  is  bounded:  find  endpoints  and,  if  X  is  below 

*  lower  endpoint  set  F(x)  =  0;  if  X  is  above  upper 

*  endpoint  set  F(x)  =  1 

XL OVER  «  LAMBDA(l)  -  1 .0/LAMBDA (2) 

IF (X . GT . XLOVER)  GOTO  10 

CDF  ■  0.0 

RETURN 

10  XUPPER  *  LAMBDA (1)  +  1 .O/LAMBDA (2) 

IF (X . LT . XUPPER)  GOTO  20 

CDF  *  1.0 

RETURN 

*  Find  the  value  of  the  CDF  at  X  using  the  secant 

*  method  to  find  p  such  that  the  inverse  CDF  evaluated 

*  at  p  is  X 

20  PI  =  0.4 
P2  =  0.6 

FI  =  GLDQF (LAMBDA, PI)  -  X 
F2  =  GLDQF (LAMBDA, P2)  -  X 
DO  30  N  -  1,  200 

PHEW  =  P2  -  (P2  -  P1)*F2/(F2  -  Fl) 

IF (PHEW. LE. 0.0)  PMEW  =  .000001 
IF(PNEW.GE.l.O)  PMEW  =  .999999 
IF(PHEW.EQ.P2)  GOTO  40 
FHEW  =  GLDqF (LAMBDA, PNEW)  -  X 
IF (DABS (FNEW).LE. 0.000005)  GOTO  40 
PI  =  P2 
P2  =  PNEW 
Fl  =  F2 

30  F2  =  FNEW 

*  Secant  method  does  not  converge 

WRITE(6,999) 

999  FORMAT (35H0ERR0R  (CDF):  SECANT  METHOD  DOES  NOT, 
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+30H  CONVERGE  AFTER  200  ITERATIONS) 

RETURN 

*  Secant  method  converges  to  a  solution 
40  CDF  -  PNEW 

RETURN 

END 

*  GLD  Quantile  Function 

FUNCTION  GLDQF  (  LAMBDA,  P  ) 

*  Declare  variables 

IMPLICIT  DOUBLE  PRECISION  (  A-H,  0-Z  ) 
DOUBLE  PRECISION  LAMBDA (4) 

*  Error  checking 

IF  (  P  .LT.  O.ODO  )  GOTO  998 
IF  (  P  .GT.  l.ODO  )  GOTO  998 

*  Compute  GLD  quantile  function 

TERM  *  P**LAMBDA(3)  -  (1.0  -  P)**LAMBDA(4) 
GLDQF  *  LAMBDA (1)  +  TERM/LAMBDA (2) 

RETURN 

*  Error  trap 

998  WRITE (6, 999)  P 

999  FORMAT (20H0ERR0R  (GLDQF):  P  *  ,  E9.4, 

+32H  AND  IS  NOT  BETWEEN  ZERO  AID  ONE) 

RETURN 

END 
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Appendix  T.  FORTRAN  Programs:  Conventional  Moment  Calculator 

***************************************************************************** 

*  FORTRAM  77  CODE  USED  II  AFIT/EMS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  CALCFX 

*  SUBROUTINE  PURPOSE:  CALCULATE  THE  OBJECTIVE  FUNCTION  FOR  POWELL’S 

*  ALGORITHM 

*  SUBROUTINE  DESCRIPTION: 

*  Evaluates  the  objective  function  to  be  minimized, 

*  which  is  the  sum  of  the  squared  differences  between 

*  the  desired  and  calculated  values  of  skewness  and 

*  kurtosis 

*  AUTHOR(S) : 

*  E.  F.  MYKYTKA 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  1979 

*  INPUT  PARAMETERS: 

*  L  -  DOUBLE  PRECISION  ARRAY  two  variables  of  function 

*  OUTPUT  PARAMETERS: 

*  FH  -  DOUBLE  PRECISION  function  value 

*  N  -  INTEGER  count  of  subroutine  calls 
***************************************************************************** 

SUBROUTINE  CALCFX (N,  L,  FN) 

*  DECLARE  LOCAL  VARIABLES  AND  GLOBAL  PARAMETERS 

IMPLICIT  REAL *8  (A-H,  0-Z) 

REAL*8  L(2) ,  LAMBDA (4) ,  MEAN 

COMMON  ALPHA3,  A3,  ALPHA4,  A4,  MEAN,  VAR,  IFUN 

*  COUNT  NUMBER  TIMES  SUBROUTINE  IS  CALLED 

IFUN  -  IFUN  +  1 

*  ENSURE  THAT  LAMBDA(3)  AND  LAMBDA (4)  HAVE  SAME  SIGN 

LAMBDA (3)  =  L(l) 

LAMBDA (4)  =  L(2) 

IF  (LAMBDA(3)*LAMBDA(4) .LT.ODO)  GOTO  999 

*  COMPUTE  OBJECTIVE  FUNCTION 

10  B40  -  l.D0/(4.D0*LAMBDA(3)+l.D0) 

B31  =  BETA (3 . DO*LAMBDA (3) +1 . DO ,  LAMBDA (4) +1 .DO) 

B22  =  BETA(2.D0*LAMBDA(3)+1.D0, 
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+  2 .D0*LAMBDA(4)+1 .DO) 

B13  -  BETA (LAMBDA (3) +1. DO,  3 .D0*LAMBDA(4)+1 .DO) 
B04  *  l.DO/(4.DO*LAMBDA(4)+l.DO) 

D  *  B40  -  4.D0*B31  +  6.DO*B22  -  4.D0*B13  +  B04 
20  B30  -  l.D0/(3.D0*LAMBDA(3)+l.D0) 

B21  «  BETA(2.D0*LAMBDA(3)+1.D0,  LAMBDA (4) +1 .DO) 
B12  *  BETA (LAMBDA (3) +1. DO,  2 . DO*LAMBDA (4) +1 . DO) 
B03  «  l.D0/(3.D0*LAMBDA(4)+l.D0) 

C  -  B30  -  3.D0*B21  +  3.D0*B12  -  B03 
30  B20  »  l.D0/(2.D0*LAMBDA(3)+l.D0) 

Bll  -  BETA (LAMBDA (3) +1 .D0,LAMBDA(4)+1 .DO) 

B02  =  1 .D0/(2 .D0*LAMBDA(4)+1 .DO) 

B  *  B20  -  2.D0*B11  +  B02 
BIO  -  1. DO/ (LAMBDA (3) +1. DO) 

B01  =  1. DO/ (LAMBDA (4) +1. DO) 

A  -  BIO  -  B01 

40  ALPHA3  ■  (C  -  3*B*A  +  2*A**3)/(B  -  A*A)**1.6 
ALPHA3  *  ALPHA3*DSIGN(1 .DO,  LAMBDA (3)) 
IF(LAMBDA(3) . NE.O.DO)  GOTO  50 
ALPHA3  =  ALPHA3*DS IGN ( 1 . DO , LAMBDA (4) ) 

50  ALPHA4  ■  D  -  4.D0*C*A  +  6.D0*B*A*A  -  3.D0*A**4 
ALPHA4  *  ALPHA4/ (B  -  A*A)**2 
DIFF3  =  ALPHA3  -  A3 
DIFF4  =  ALPHA4  -  A4 

*  FUNCTION  VALUE  IS  SQUARED  DIFFERENCE  BETWEEN 

*  CALCULATED  AND  DESIRED  ALPHA3  AND  ALPHA4 

FN  *  DIFF3*DIFF3  +  DIFF4*DIFF4 
RETURN 

999  FN  =  10. DO 
RETURN 
END 


♦*****♦**♦♦ ********* ************************** **♦♦****♦**♦♦***♦* ************* 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  CMOMENTS 

*  SUBROUTINE  PURPOSE:  CALCULATE  THE  FIRST  FOUR  CONVENTIONAL  MOMENTS  OF  A 

*  DATA  SAMPLE 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  receives  an  array  containing  the  data  sample,  and  computes 

*  the  mean,  variance,  skewness,  and  kurtosis  of  that  sample.  It  then 
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*  returns  those  values  in  an  output  parameter  array. 

*  AUTHOR(S) : 

*  R.  B.  HORAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  VRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  NOVEMBER  1993 

*  INPUT  PARAMETERS: 

*  ARRAY  -  DOUBLE  PRECISION  ARRAY  sample  data 

*  LENGTH  -  INTEGER  length  of  input  array  (number  of  elements  in  sample) 

*  OUTPUT  PARAMETERS: 

*  MOMENT  -  DOUBLE  PRECISION  ARRAY  first  four  conventional  moments 
***************************************************************************** 

SUBROUTINE  CMOMENTS  (  ARRAY,  LENGTH,  MOMENT  ) 

*  DECLARE  VARIABLES 

INTEGER  LENGTH,  I 

DOUBLE  PRECISION  ARRAY  (  1:3000  ) 

DOUBLE  PRECISION  MOMENT  (  1:4  ) 

DOUBLE  PRECISION  MEAN,  VAR,  SKEW,  KURT 

*  INITIALIZE  VARIABLES 

MEAN  =  0 
VAR  =  0 
SKEW  =  0 
KURT  ■  0 

*  COMPUTE  MEAN 

DO  101  I  =  1,  LENGTH 
MEAN  =  MEAN  +  ARRAY  (I) 

101  CONTINUE 

MEAN  =  MEAN/LENGTH 

*  COMPUTE  VARIANCE,  SKEWNESS,  AND  KURTOSIS  SUMS 

DO  102  1=1,  LENGTH 

VAR  =  VAR  +  (ARRAY (I) -MEAN) **2 
SKEW  =  SKEW  +  (ARRAY(I)-MEAN)**3 
KURT  =  KURT  +  (ARRAY (I) -MEAN) **4 

102  CONTINUE 

*  DIVIDE  SUMS  BY  NUMBER  OF  ELEMENTS  IN  SAMPLE 

VAR  =  VAR/LENGTH 
SKEW  =  SKEW/LENGTH 
KURT  =  KURT/LENGTH 

*  COMPIT  Z  THIRD  AND  FOURTH  MOMENTS 

SKEW  =  SKEW/ (VAR* *1.5) 
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KURT  *  KURT/ (VAR**2) 

*  LOAD  OUTPUT  PARAMETER  ARRAY 

MOMENT(l)  =  MEAN 
MOMENT (2)  =  VAR 
MOMENT (3)  =  SKEW 
MOMENT (4)  =  KURT 

*  RETURN  TO  PROGRAM 

RETURN 

END 
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Appendix  U.  FORTRAN  Programs:  Alternate  Moment  Calculator 

***************************************************************************** 

*  FORTRAH  77  CODE  USED  IH  AFIT/EHS  MASTER'S  THESIS; 

*  "THE  USE  OF  L-HOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  CALCFX 

*  SUBROUTINE  PURPOSE:  CALCULATE  THE  OBJECTIVE  FUNCTION  FOR  POWELL'S 

*  ALGORITHM 

*  SUBROUTINE  DESCRIPTION: 

*  Evaluates  the  objective  function  to  be  minimized, 

*  which  is  the  sum  of  the  squared  differences  between 

*  the  desired  and  calculated  values  of  skewness  and 

*  kurtosis 

*  AUTHOR (S) : 

*  E.  F.  MYKYTKA 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  1979 

*  INPUT  PARAMETERS: 

*  L  -  DOUBLE  PRECISION  ARRAY  two  variables  of  function 

*  OUTPUT  PARAMETERS: 

*  FN  -  DOUBLE  PRECISION  function  value 

*  N  -  INTEGER  count  of  subroutine  calls 

SUBROUTINE  CALCFX  (N,  L,  FN) 

*  DECLARE  LOCAL  VARIABLES  AND  GLOBAL  PARAMETERS 

IMPLICIT  REAL*8  (A-H,  0-Z) 

REAL*8  L(2)  ,  LAMBDA (4)  ,  MEAN 

COMMON  ALPHA3,  A3,  ALPHA4,  A4,  MEAN,  VAR,  I FUN 

*  COUNT  NUMBER  TIMES  SUBROUTINE  IS  CALLED 

I FUN  *  IFUN  +  1 

*  ENSURE  THAT  LAMBDA (3)  AND  LAMBDA (4)  HAVE  SAME  SIGN 

L3  ■  L(1)+1.0D0 
L4  =  L(2)+1.0D0 

IF  ((L(1)*L(2))  .LT.  O.ODO)  GOTO  999 

*  COMPUTE  OBJECTIVE  FUNCTION 

U06-20 . DO* ( ( (1 .DO- (0 . 95D0**L3) ) /L3) - ( (0 . 05D0**L4) /L4) ) 

U60=2.D0*  ( ( (1 . DO- (0 . 5D0**L3  ))/L3)-( (0.5D0**L4  )/L4)) 

L05*»20 .  DO*  (  (  (0 . 05D0**L3)  /L3) + (  (  (0 . 95D0**L4)  -1 .  DO)  /L4)  ) 
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L50-2.D0*  (((0.5D0**L3  )/L3)+(((0.5D0**L4  )~1.D0)/L4)) 

M50= ( (0 . 76D0**L3) - (0 . 26D0**L3) ) /L3 

M50=2 . DO* (M60+ ( ( (0 . 25D0**L4) - (0 . 75D0**L4) ) /L4) ) 

ALPHA3= (U05-M50) / (M60-L05) 

ALPHA4= (U05-L05) / (USO-LSO) 

DIFF3  =  ALPHAS  -  A3 
DIFF4  *  ALPHA4  -  A4 

*  FUNCTION  VALUE  IS  SQUARED  DIFFERENCE  BETWEEN 

*  CALCULATED  AND  DESIRED  ALPHA3  AND  ALPHA4 

FN  -  DIFF3*DIFF3  +  DIFF4*DIFF4 
RETURN 

999  FH  =  10. DO 
RETURN 
END 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER'S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  AMOHENTS 

*  SUBROUTINE  PURPOSE:  CALCULATE  THE  FIRST  FOUR  ALTERNATE  MOMENTS  OF  A 

*  DATA  SAMPLE 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  receives  an  array  containing  the  sorted  data  sample ,  and 

*  computes  the  mean,  variance,  Q3,  and  Q4  of  that  sample.  It  then 

*  returns  those  values  in  an  output  parameter  array. 

*  AUTHOR(S) : 

*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  NOVEMBER  1993 

*  INPUT  PARAMETERS: 

*  ARRAY  -  DOUBLE  PRECISION  ARRAY  sample  data 

*  LENGTH  -  INTEGER  length  of  input  array  (number  of  elements  in  sample) 

*  OUTPUT  PARAMETERS : 

*  MOMENT  -  DOUBLE  PRECISION  ARRAY  first  four  alternate  moments 
***************************************************************************** 

SUBROUTINE  AUGMENTS  (  ARRAY,  LENGTH,  MOMENT  ) 

*  DECLARE  VARIABLES 

INTEGER  LENGTH,  I,  A,  B 
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DOUBLE  PRECISION  ARRAY  (  1:3000  ) 

DOUBLE  PRECISION  MOMENT  (  1:4  ) 

DOUBLE  PRECISION  MEAN,  VAR,  U,  L,  M,  UH,  LH 

*  INITIALIZE  VARIABLES 

MEAN  »  0 
VAR  *  0 
A  =  1 
B  -  12 
U  -  0 
L  =  0 
M  -  0 
UH  -  0 
LH  =  0 

*  CALCULATE  MEAN 

DO  101  1=1,  LENGTH 
MEAN  =  MEAN  +  ARRAY  (I) 

101  CONTINUE 

MEAN  =  MEAN/LENGTH 

*  CALCULATE  VARIANCE 

DO  102  1*1,  LENGTH 
VAR  *  VAR  +  (ARRAY (I) -MEAN) **2 

*  CALCULATE  q3  AND  Q4  (SAMPLE  SIZE  OF  25  ONLY) 

102  CONTINUE 

U  =  U  +  ARRAY (25)  +  (ARRAY(24)/4) 

L  -  L  +  ARRAY(l)  +  (ARRAY(2)/4) 

DO  104  I  -  1,  B 
UH  *  UH  +  ARRAY (1+13) 

LH  -  LH  +  ARRAY (I) 

104  CONTINUE 

UH  ■  UH  +  (ARRAY(13)/2) 

LH  -  LH  +  (ARRAY(13)/2) 

DO  105  I  -  8,  18 
M  =  M  +  ARRAY(I) 

105  CONTINUE 

M  ■  M  +  (ARRAY (7)  *  0.75D0) 

M  =  M  +  (ARRAY(19)  *  0.75D0) 

U=U/1 . 25 
L-L/1.25 
UH=UH/12 . 5 
LH-LH/12.5 
M=M/12.5 

*  LOAD  OUTPUT  PARAMETER  ARRAY 
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NONEHT(l)  -  MEAN 
MOMENT (2)  »  VAR/LENGTH 
M0MEMT(3)  »  (U-M) / (M-L) 

MOMENT (4)  »  (U-L) / (UH-LH) 

RETURN  TO  PROGRAM 
RETURN 
END 
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Appendix  V.  FORTRAN  Program:  Output  File  Statistic  Calculator 

***************************************************************************** 

*  FORTRAN  77  CODE  USED  IN  AFIT/ENS  MASTER’S  THESIS; 

*  "THE  USE  OF  L-MOMENTS  TO  FIT  THE  GENERALIZED  LAMBDA  DISTRIBUTION 

*  TO  SAMPLE  DATA" 

*  SUBROUTINE  TITLE:  STAT 

*  SUBROUTINE  PURPOSE: 

*  SUBROUTINE  DESCRIPTION: 

*  This  subroutine  receives  the  output  from  30  reports  and  finds  the  minimum 

*  and  maximum  values  of  each  of  11  data  types  in  the  reports.  It  then 

*  calculates  the  mean  and  variance  of  the  30  values  for  each  data  type. 

*  Then  it  outputs  the  statistics  in  a  LaTex  tabular  environment  format. 

*  AUTHOR(S) : 

*  R.  B.  MOHAN,  CAPTAIN,  USAF 

*  DEPARTMENT  OF  OPERATIONAL  SCIENCES 

*  GRADUATE  SCHOOL  OF  ENGINEERING,  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

*  WRIGHT-PATTERSON  AIR  FORCE  BASE 

*  OHIO  45433,  USA 

*  DATE:  NOVEMBER  1993 

*  INPUT  PARAMETERS: 

*  User  inputs  name  of  input  data  file 

*  OUTPUT  PARAMETERS: 

*  User  inputs  name  of  output  statistics  file 
***************************************************************************** 

PROGRAM  STAT 

******  DECLARE  LOCAL  VARIABLES 

DOUBLE  PRECISION  DATA  (1:11,  1:30) 

*  11  values  by  30  samples 

DOUBLE  PRECISION  STATS  (1:11,  1:4) 

*  11  values  by  4  statistics 

INTEGER  I,  J 

DOUBLE  PRECISION  MIN,  MAX,  AVG,  VAR 
CHARACTER* 12  INFILE,  OUTFILE 
******  PROMPT  USER  FOR  NAME  OF  FILE 

PRINT  *,  ’ENTER  THE  NAME  OF  THE  INPUT  FILE:’ 

READ  (*,  ’ (A12) ’  )  INFILE 
******  OPEN  INPUT  FILE  AND  UPLOAD  FORMATTED  DATA 
OPEN  (UNIT=20 ,  FILE=INFILE,  STATUS*’ OLD’ ) 

10  FORMAT  (F13.6) 

******  download  input  file  to  array 
DO  100  I  =  1,  30 
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DO  120  J  -  1,  11 

READ  (20,  10,  END=110  )  DATA  (J,I) 

120  CONTINUE 
ioo  com  hue 

******  close  output  file 
110  CLOSE  (20) 

******  CALCULATE  STATISTICS 
DO  200  1*1,  11 

*  INITIALIZE  VARIABLES 

MIH  *  10000 
MAI  -  -10000 
AVG  -  0 
VAR  =  0 

*  FIND  MIN,  MAX,  AND  SUM  OF  VALUES 

DO  300  J  -  1,  30 

IF  (DATA (I, J)  .LT.  MIN)  MIN  =  DATA(I,J) 

IF  (DATA(I.J)  .GT.  MAI)  MAX  =  DATA(I.J) 

AVG  -  AVG  +  DATA(I,J) 

300  CONTINUE 

*  DIVIDE  SUM  BY  NUMBER  OF  VALUES  TO  FIND  MEAN 

AVG  =  AVG/30 

*  CALCULATE  VARIANCE 

DO  400  J  =  1,  30 

VAR  *  VAR  +  (AVG-DATA(I,J))**2 
400  CONTINUE 

*  LOAD  OUTPUT  PARAMETER  ARRAY 

STATS (1,1)  =  MIN 
STATS (I, 2)  *  MAX 
STATS (I, 3)  *  AVG 
STATS (I, 4)  =  VAR/29 
200  CONTINUE 

******  PROMPT  USER  FOR  NAME  OF  FILE 

PRINT  *,  ’ENTER  THE  NAME  OF  THE  OUTPUT  FILE:’ 

READ  (*,  ’ (A12) ’  )  OUTFILE 
******  open  OUTPUT  FILE  AND  DOWNLOAD  FORMATTED  DATA 
OPEN  (UNIT=30 ,  FILE=OUTFILE,  STATUS* ’ NEW ’ ) 

******  FORMAT  FOR  LATEX  TABULAR  ENVIRONMENT 

40  FORMAT  (’  ft  1st  ft  \F11.8,’  ft  ’, 

+F11.8,’  ft  ’ ,F11 .8, ’  ft  ’ ,F11 .8, ’  \\\\  \\cline{3-6}’) 

41  FORMAT  (’Moments  of  ft  2nd  ft  ’,F11.8,’  ft  ', 

+F11.8,’  ft  ’ ,F11 .8, ’  ft  ’ ,F11 .8, ’  \\\\  \\cline{3-6}’) 

42  FORMAT  (’Sample  Data  ft  3rd  ft  ’.F11.8,’  ft 
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+F11.8,’  ft  ’.F11.8,’  ft  ’.F11.8,’  \\\\  \\cline{3-6}’) 

43  FORMAT  (’  ft  4th  ft  \Fli.8,’  ft  ’, 

+F11 .8,  *  ft  \F11.8,’  ft  *  ,F11 .8,  *  \\\\  Whline’) 

44  FORMAT  (’  ft  \\1  ft  1 ,F11 .8, ’  ft  ’ , 

+F11.8,’  »  \F11.8,’  ft  ’ ,F11 .8, ’  \\\\  \\cline{3-6}’) 

45  FORMAT  (’Lambda  ft  \\2  ft  ’.F11.8,’  ft  ’, 

+F11 .8, ’  ft  »,F11.8,»  ft  \F11.8,’  \\\\  \\cline{3-6}’) 

46  FORMAT  (’Parameters  ft  \\3  ft  »,F11.8,’  ft  ’, 

+F11.8,*  ft  »,FU.8,»  ft  *  ,F11 .8,  *  \\\\  \\cline{3-6}>) 

47  FORMAT  (»  ft  \\4  ft  » ,F11.8,»  ft  \ 

+F11.8,’  ft  * ,F11 .8, *  ft  »,F11.8,»  \\\\  Whline’) 

48  FORMAT  (’Kolmogorov-  ft  MIN  ft  \F11.8,’  ft  ’, 

+F11.8,’  ft  ’ ,F11 .8, ’  ft  ’ ,F11 .8, ’  \\\\  \\cline{3-6>’) 

49  FORMAT  (’Smirnov  ft  AVG  ft  ’.F11.8,’  ft  ’, 

+F11.8,’  ft  ' ,F11 .8, '  ft  ’ ,F11 .8, ’  \\\\  \\cline{3-6}’) 

50  FORMAT  (’Statistics  ft  MAX  ft  ’.F11.8,’  ft  ’, 

+F11.8,’  ft  ’  ,F11 .8,  ’  ft  ’  ,F11 .8,  ’  \\\\  WhlineWhline’) 

******  download  array  to  output  file 

HRITE  (30,40)  STATS(l.l),  STATS(l,2),  STATS(1,3),  STATS(1,4) 
WRITE  (30,41)  STATS (2,1) ,  STATS(2,2),  STATS(2,3),  STATS(2,4) 
WRITE  (30,42)  STATS (3,1),  STATS (3, 2),  STATS (3, 3),  STATS (3, 4) 
WRITE  (30,43)  STATS (4,1),  STATS (4, 2),  STATS (4, 3),  STATS (4, 4) 
WRITE  (30,44)  STATS(5,1),  STATS(5,2),  STATS(5,3),  STATS(5,4) 
WRITE  (30,45)  STATS (6,1),  STATS (6, 2),  STATS (6, 3),  STATS (6, 4) 
WRITE  (30,46)  STATS(7,1),  STATS(7,2),  STATS(7,3),  STATS(7,4) 
WRITE  (30,47)  STATS(8,1),  STATS(8,2),  STATS(8,3),  STATS(8,4) 
WRITE  (30,48)  STATS (9,1),  STATS (9, 2),  STATS (9, 3),  STATS (9, 4) 
WRITE  (30,49)  STATS(10,1),  STATS(10,2),  STATS(10,3),  STATS(10,4) 
WRITE  (30,50)  STATS (11 ,1) ,  STATS(11,2),  STATS(11,3),  STATS(11,4) 
******  close  output  file 
CLOSE  (30) 

STOP 

END 
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The  Generalized  Lambda  Distribution  (GLD)  is  a  four-parameter,  continuous  probability  distribution  that  is 
useful  for  simulation  analysis.  The  strengths  of  the  GLD  lie  in  its  abilities  to  approximate  many  distributions, 
represent  data  when  the  underlying  distribution  is  unknown,  and  fit  or  generate  random  variates.  The  method 
of  momenta  is  presently  the  accepted  technique  for  estimating  the  parameters  of  this  distribution.  However,  ^ 

it  is  sensitive  to  extreme  observations  and  subject  to  large  sampling  variability  as  the  sample  size  decreases. 
L-moments  are  expectations  of  certain  linear  combinations  of  order  statistics.  They  can  be  used  to  estimate 
parameters  and  quantiles  of  probability  distributions.  Their  main  advantage  over  conventional  moments  is 
that  they  suffer  less  from  the  effects  of  sampling  variability,  and  are  theoretically  more  robust  to  outliers 
than  conventional  moments.  Estimating  the  parameters  of  the  GLD  by  matching  its  L-moments  to  those 
of  the  sample  is  known  as  the  method  of  L-moments.  This  appears  to  be  an  attravtive  alternative  to  the 
method  of  moments  and  is  developed  in  this  thesis.  A  Monte  Carlo  experiment  compared  the  method  of 
L-moments  to  the  methed  of  conventional  moments  and  a  third  method  which  uses  alternate  measures  of 
symmetry  and  tailweight.  Experiment  results  showed  that  L-moments  are  better  than  conventional  and 
alternate  moments  for  fitting  distributions  to  sample  data,  particularly  when  the  skewness  and  kurtosis  of 
the  sample  distribution  are  large. 
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